A scalable parallel algorithm for large-scale reactive force-field molecular dynamics simulations

A scalable parallel algorithm has been designed to perform multimillion-atom molecular dynamics (MD) simulations, in which first principles-based reactive force fields (ReaxFF) describe chemical reactions. Environment-dependent bond orders associated with atomic pairs and their derivatives are reused extensively with the aid of linked-list cells to minimize the computation associated with atomic n-tuple interactions (n?4 explicitly and ?6 due to chain-rule differentiation). These n-tuple computations are made modular, so that they can be reconfigured effectively with a multiple time-step integrator to further reduce the computation time. Atomic charges are updated dynamically with an electronegativity equalization method, by iteratively minimizing the electrostatic energy with the charge-neutrality constraint. The ReaxFF-MD simulation algorithm has been implemented on parallel computers based on a spatial decomposition scheme combined with distributed n-tuple data structures. The measured parallel efficiency of the parallel ReaxFF-MD algorithm is 0.998 on 131,072 IBM BlueGene/L processors for a 1.01 billion-atom RDX system.     

Stochastic kinetic Monte Carlo algorithms for long-range Hamiltonians

We present a higher order kinetic Monte Carlo methodology suitable to model the evolution of systems in which the transition rates are non-trivial to calculate or in which Monte Carlo moves are likely to be non-productive flicker events. The second order residence time algorithm first introduced by Athènes et al. [Phil. Mag. A 76 (1997) 565] is rederived from the n-fold way algorithm of Bortz et al. [J. Comput. Phys. 17 (1975) 10] as a fully stochastic algorithm. The second order algorithm can be dynamically called when necessary to eliminate unproductive flickering between a metastable state and its neighbours. An algorithm combining elements of the first order and second order methods is shown to be more efficient, in terms of the number of rate calculations, than the first order or second order methods alone while remaining statistically identical. This efficiency is of prime importance when dealing with computationally expensive rate functions such as those arising from long-range Hamiltonians. Our algorithm has been developed for use when considering simulations of vacancy diffusion under the influence of elastic stress fields. We demonstrate the improved efficiency of the method over that of the n-fold way in simulations of vacancy diffusion in alloys. Our algorithm is seen to be an order of magnitude more efficient than the n-fold way in these simulations. We show that when magnesium is added to an Al-2at.%Cu alloy, this has the effect of trapping vacancies. When trapping occurs, we see that our algorithm performs thousands of events for each rate calculation performed.     

xPerm: fast index canonicalization for tensor computer algebra

We present a very fast implementation of the Butler-Portugal algorithm for index canonicalization with respect to permutation symmetries. It is called xPerm, and has been written as a combination of a Mathematica package and a C subroutine. The latter performs the most demanding parts of the computations and can be linked from any other program or computer algebra system. We demonstrate with tests and timings the effectively polynomial performance of the Butler-Portugal algorithm with respect to the number of indices, though we also show a case in which it is exponential. Our implementation handles generic tensorial expressions with several dozen indices in hundredths of a second, or one hundred indices in a few seconds, clearly outperforming all other current canonicalizers. The code has been already under intensive testing for several years and has been essential in recent investigations in large-scale tensor computer algebra.

Program summary

Program title: xPermCatalogue identifier: AEBH_v1_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/AEBH_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.: 93 582No. of bytes in distributed program, including test data, etc.: 1 537 832Distribution format: tar.gzProgramming language: C and Mathematica (version 5.0 or higher)Computer: Any computer running C and Mathematica (version 5.0 or higher)Operating system: Linux, Unix, Windows XP, MacOSRAM:: 20 MbyteWord size: 64 or 32 bitsClassification: 1.5, 5Nature of problem: Canonicalization of indexed expressions with respect to permutation symmetries.Solution method: The Butler-Portugal algorithm.Restrictions: Multiterm symmetries are not considered.Running time: A few seconds with generic expressions of up to 100 indices. The xPermDoc.nb notebook supplied with the distribution takes approximately one and a half hours to execute in full.     

Scalable and portable visualization of large atomistic datasets

A scalable and portable code named Atomsviewer has been developed to interactively visualize a large atomistic dataset consisting of up to a billion atoms. The code uses a hierarchical view frustum-culling algorithm based on the octree data structure to efficiently remove atoms outside of the user's field-of-view. Probabilistic and depth-based occlusion-culling algorithms then select atoms, which have a high probability of being visible. Finally a multiresolution algorithm is used to render the selected subset of visible atoms at varying levels of detail. Atomsviewer is written in C++ and OpenGL, and it has been tested on a number of architectures including Windows, Macintosh, and SGI. Atomsviewer has been used to visualize tens of millions of atoms on a standard desktop computer and, in its parallel version, up to a billion atoms.

Program summary

Title of program: AtomsviewerCatalogue identifier: ADUMProgram summary URL:http://cpc.cs.qub.ac.uk/summaries/ADUMProgram obtainable from: CPC Program Library, Queen's University of Belfast, N. IrelandComputer for which the program is designed and others on which it has been tested: 2.4 GHz Pentium 4/Xeon processor, professional graphics card; Apple G4 (867 MHz)/G5, professional graphics cardOperating systems under which the program has been tested: Windows 2000/XP, Mac OS 10.2/10.3, SGI IRIX 6.5Programming languages used: C++, C and OpenGLMemory required to execute with typical data: 1 gigabyte of RAMHigh speed storage required: 60 gigabytesNo. of lines in the distributed program including test data, etc.: 550 241No. of bytes in the distributed program including test data, etc.: 6 258 245Number of bits in a word: ArbitraryNumber of processors used: 1Has the code been vectorized or parallelized: NoDistribution format: tar gzip fileNature of physical problem: Scientific visualization of atomic systemsMethod of solution: Rendering of atoms using computer graphic techniques, culling algorithms for data minimization, and levels-of-detail for minimal renderingRestrictions on the complexity of the problem: NoneTypical running time: The program is interactive in its executionUnusual features of the program: NoneReferences: The conceptual foundation and subsequent implementation of the algorithms are found in [A. Sharma, A. Nakano, R.K. Kalia, P. Vashishta, S. Kodiyalam, P. Miller, W. Zhao, X.L. Liu, T.J. Campbell, A. Haas, Presence—Teleoperators and Virtual Environments 12 (1) (2003)].     

Fast clear-sky solar irradiation computation for very large digital elevation models

This paper presents a fast algorithm to compute the global clear-sky irradiation, appropriate for extended high-resolution Digital Elevation Models (DEMs). The latest equations published in the European Solar Radiation Atlas (ESRA) have been used as a starting point for the proposed model and solved using a numerical method. A new calculation reordering has been performed to (1) substantially diminish the computational requirements, and (2) to reduce dependence on both, the DEM size and the simulated period, i.e., the period during which the irradiation is calculated. All relevant parameters related to shadowing, atmospheric, and climatological factors have been considered. The computational results demonstrate that the obtained implementation is faster by many orders of magnitude than all existing advanced irradiation models while maintaining accuracy. Although this paper focuses on the clear-sky irradiation, the developed software also computes the global irradiation applying a filter that considers the clear-sky index.     

dHybrid: A massively parallel code for hybrid simulations of space plasmas

A massively parallel simulation code, called dHybrid, has been developed to perform global scale studies of space plasma interactions. This code is based on an explicit hybrid model; the numerical stability and parallel scalability of the code are studied. A stabilization method for the explicit algorithm, for regions of near zero density, is proposed. Three-dimensional hybrid simulations of the interaction of the solar wind with unmagnetized artificial objects are presented, with a focus on the expansion of a plasma cloud into the solar wind, which creates a diamagnetic cavity and drives the Interplanetary Magnetic Field out of the expansion region. The dynamics of this system can provide insights into other similar scenarios, such as the interaction of the solar wind with unmagnetized planets.     

A new algorithm for fast simulation of optical spectra generated in a nonlinear media

A new and highly efficient algorithm developed under MATLAB for calculating the optical spectra generated by non-resonant optical parametric fluorescence is presented. This algorithm, which allows quick simulation of the spectra, is shown to be much more rapid than standard ones. The ways to modify the algorithm for other environments are discussed.

Program summary

Title of the program:OPFSpectraFinderCatalogue number: ADTKProgram summary URL:http://cpc.cs.qub.ac.uk/summaries/ADTKProgram obtainable from: CPC Program Library, Queen's University of Belfast, N. IrelandLicensing provisions: Persons requesting the program must sign the standard CPC non-profit use licenseComputers on which the program has been tested: Intel Pentium, SUN HPC3500Operating system under which the program, has been tested: Microsoft Windows 2000, Unix Solaris 7.0Programming language used: MATLAB 6.0Memory required to execute with typical data: 256 M (the total computer memory)No. of processors used: 1No. of bytes in distributed program, including test data, etc.: 935 376No. of lines in distributed program, including test data, etc.: 6667Distribution format: tar gzip fileNature of physical problem: Obtaining the form of optical spectra of light generated in the process of optical parametric fluorescence in a periodically poled quadratically nonlinear mediumMethod of solution: Linearization of the set of nonlinear equations, obtaining the set of equations with constant coefficients, finding the solution as the matrix exponentialRestrictions on the complexity of the problem: not knownTypical running time: 80 seconds for obtaining the data for data set by default in the files submitted on a PC mentioned aboveReferences: The physical discussion is given in V. Beskrovnyy, P. Baldi, Optical parametric fluorescence spectra in periodically poled media, Optics Express 10 (2002) 990     

Convergence and refinement of the Wang-Landau algorithm

Recently, Wang and Landau proposed a new random walk algorithm that can be very efficiently applied to many problems. Subsequently, there has been numerous studies on the algorithm itself and many proposals for improvements were put forward. However, fundamental questions such as what determines the rate of convergence has not been answered. To understand the mechanism behind the Wang-Landau method, we did an error analysis and found that a steady state is reached where the fluctuations in the accumulated energy histogram saturate at values proportional to [log(f)]^−1/2. This value is closely related to the error corrections to the Wang-Landau method. We also study the rate of convergence using different “tuning” parameters in the algorithm.     

Large scale fractal aggregates using the tunable dimension cluster-cluster aggregation

The tunable dimension cluster-cluster aggregation (tdCCA) [R. Thouy, R. Jullien, J. Phys. A: Math. Gen. 27 (1994) 2953] provides a computational model for creating fractal aggregates with a tunable fractal dimension. A straightforward implementation of this model requires a computational effort scaling with O(N_total⁴) of the number of particles N_total. By applying two minor changes to the algorithm the computational effort can be reduced to O(N_total²) and allows an efficient parallel implementation of the tdCCA. On a modern parallel computer a fractal aggregate of one million particles has been built in less than 24 h.     

Symbolic computation of hyperbolic tangent solutions for nonlinear differential-difference equations

A new algorithm is presented to find exact traveling wave solutions of differential-difference equations in terms of tanh functions. For systems with parameters, the algorithm determines the conditions on the parameters so that the equations might admit polynomial solutions in tanh. Examples illustrate the key steps of the algorithm. Through discussion and example, parallels are drawn to the tanh-method for partial differential equations. The new algorithm is implemented in Mathematica. The package DDESpecialSolutions.m can be used to automatically compute traveling wave solutions of nonlinear polynomial differential-difference equations. Use of the package, implementation issues, scope, and limitations of the software are addressed.

Program summary

Title of program: DDESpecialSolutions.mCatalogue identifier:ADUJProgram summary URL:http://cpc.cs.qub.ac.uk/summaries/ADUJProgram obtainable from:CPC Program Library, Queen's University of Belfast, N. IrelandDistribution format: tar.gzComputers: Created using a PC, but can be run on UNIX and Apple machinesOperating systems under which the program has been tested: Windows 2000 and Windows XPProgramming language used: Mathematica, version 3.0 or higherMemory required to execute with typical data: 9 MBNumber of processors used: 1Has the code been vectorised or parallelized?: NoNumber of lines in distributed program, including test data, etc.: 3221Number of bytes in distributed program, including test data, etc.: 23 745Nature of physical problem: The program computes exact solutions to differential-difference equations in terms of the tanh function. Such solutions describe particle vibrations in lattices, currents in electrical networks, pulses in biological chains, etc.Method of solution: After the differential-difference equation is put in a traveling frame of reference, the coefficients of a candidate polynomial solution in tanh are solved for. The resulting traveling wave solutions are tested by substitution into the original differential-difference equation.Restrictions on the complexity of the program: The system of differential-difference equations must be polynomial. Solutions are polynomial in tanh.Typical running time: The average run time of 16 cases (including the Toda, Volterra, and Ablowitz-Ladik lattices) is 0.228 seconds with a standard deviation of 0.165 seconds on a 2.4 GHz Pentium 4 with 512 MB RAM running Mathematica 4.1. The running time may vary considerably, depending on the complexity of the problem.     

On servers in series with losses in descrete time

We study a m-phase queueing system without buffers, operating in discrete time. The input flow is Bernoulli with parameter a. Service times in server i have geometric distribution with parameter b _i. A customer, trying to enter a server at an instant, when it is busy, is lost. There have been obtained system of equilibrium equations and recurrence relations for its coefficients which enable us to formulate the algorithm to build the system. Recurrence formulas for computation of the empty system probability and some other performance characteristics of the system, are determined. The problem of optimal allocation of the servers is studied numerically.     

An events based algorithm for distributing concurrent tasks on multi-core architectures

In this paper, a programming model is presented which enables scalable parallel performance on multi-core shared memory architectures. The model has been developed for application to a wide range of numerical simulation problems. Such problems involve time stepping or iteration algorithms where synchronization of multiple threads of execution is required. It is shown that traditional approaches to parallelism including message passing and scatter-gather can be improved upon in terms of speed-up and memory management. Using spatial decomposition to create orthogonal computational tasks, a new task management algorithm called H-Dispatch is developed. This algorithm makes efficient use of memory resources by limiting the need for garbage collection and takes optimal advantage of multiple cores by employing a “hungry” pull strategy. The technique is demonstrated on a simple finite difference solver and results are compared to traditional MPI and scatter-gather approaches. The H-Dispatch approach achieves near linear speed-up with results for efficiency of 85% on a 24-core machine. It is noted that the H-Dispatch algorithm is quite general and can be applied to a wide class of computational tasks on heterogeneous architectures involving multi-core and GPGPU hardware.     

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.     

New version of PLNoise: a package for exact numerical simulation of power-law noises

In a recent paper I have introduced a package for the exact simulation of power-law noises and other colored noises [E. Milotti, Comput. Phys. Comm. 175 (2006) 212]: in particular, the algorithm generates 1/fα noises with 0<α?2. Here I extend the algorithm to generate 1/fα noises with 2<α?4 (black noises). The method is exact in the sense that it produces a sampled process with a theoretically guaranteed range-limited power-law spectrum for any arbitrary sequence of sampling intervals, i.e. the sampling times may be unevenly spaced.

Program summary

Title of program: PLNoiseCatalogue identifier:ADXV_v2_0Program summary URL:http://cpc.cs.qub.ac.uk/summaries/ADXV_v2_0.htmlLicensing provisions: Standard CPC licence, http://cpc.cs.qub.ac.uk/licence/licence.htmlProgram obtainable from: CPC Program Library, Queen's University of Belfast, N. IrelandProgramming language used: ANSI CComputer: Any computer with an ANSI C compiler: the package has been tested with gcc version 3.2.3 on Red Hat Linux 3.2.3-52 and gcc version 4.0.0 and 4.0.1 on Apple Mac OS X-10.4Operating system: All operating systems capable of running an ANSI C compilerRAM: The code of the test program is very compact (about 60 Kbytes), but the program works with list management and allocates memory dynamically; in a typical run with average list length 2⋅¹⁰⁴, the RAM taken by the list is 200 KbytesExternal routines: The package needs external routines to generate uniform and exponential deviates. The implementation described here uses the random number generation library ranlib freely available from Netlib [B.W. Brown, J. Lovato, K. Russell: ranlib, available from Netlib, http://www.netlib.org/random/index.html, select the C version ranlib.c], but it has also been successfully tested with the random number routines in Numerical Recipes [W.H. Press, S.A. Teulkolsky, W.T. Vetterling, B.P. Flannery, Numerical Recipes in C: The Art of Scientific Computing, second ed., Cambridge Univ. Press., Cambridge, 1992, pp. 274-290]. Notice that ranlib requires a pair of routines from the linear algebra package LINPACK, and that the distribution of ranlib includes the C source of these routines, in case LINPACK is not installed on the target machine.No. of lines in distributed program, including test data, etc.:2975No. of bytes in distributed program, including test data, etc.:194 588Distribution format:tar.gzCatalogue identifier of previous version: ADXV_v1_0Journal reference of previous version: Comput. Phys. Comm. 175 (2006) 212Does the new version supersede the previous version?: YesNature of problem: Exact generation of different types of colored noise.Solution method: Random superposition of relaxation processes [E. Milotti, Phys. Rev. E 72 (2005) 056701], possibly followed by an integration step to produce noise with spectral index >2.Reasons for the new version: Extension to 1/fα noises with spectral index 2<α?4: the new version generates both noises with spectral with spectral index 0<α?2 and with 2<α?4.Summary of revisions: Although the overall structure remains the same, one routine has been added and several changes have been made throughout the code to include the new integration step.Unusual features: The algorithm is theoretically guaranteed to be exact, and unlike all other existing generators it can generate samples with uneven spacing.Additional comments: The program requires an initialization step; for some parameter sets this may become rather heavy.Running time: Running time varies widely with different input parameters, however in a test run like the one in Section 3 in the long write-up, the generation routine took on average about 75 μs for each sample.     

Kinematical calculations of RHEED intensity oscillations during the growth of thin epitaxial films

A practical computing algorithm working in real time has been developed for calculating the reflection high-energy electron diffraction (RHEED) from the molecular beam epitaxy (MBE) growing surface. The calculations are based on the use of kinematical diffraction theory. Simple mathematical models are used for the growth simulation in order to investigate the fundamental behaviors of reflectivity change during the growth of thin epitaxial films prepared using MBE.

Program summary

Title of program:GROWTHCatalogue identifier:ADVLProgram summary URL:http://cpc.cs.qub.ac.uk/summaries/ADVLProgram obtainable from: CPC Program Library, Queen's University of Belfast, N. IrelandDistribution format: tar.gzComputer for which the program is designed and others on which is has been tested:Pentium-based PCOperating systems or monitors under which the program has been tested:Windows 9x, XP, NTProgramming language used:Object PascalMemory required to execute with typical data:more than 1 MBNumber of bits in a word: 64 bitsNumber of processors used: 1Number of lines in distributed program, including test data, etc.: 10 989Number of bytes in distributed program, including test data, etc.:103 048Nature of the physical problem:Reflection high-energy electron diffraction (RHEED) is a very useful technique for studying growth and surface analysis of thin epitaxial structures prepared using the molecular beam epitaxy (MBE). The simplest approach to calculating the RHEED intensity during the growth of thin epitaxial films is the kinematical diffraction theory (often called kinematical approximation), in which only a single scattering event is taken into account. The biggest advantage of this approach is that we can calculate RHEED intensity in real time. Also, the approach facilitates intuitive understanding of the growth mechanism and surface morphology [P.I. Cohen, G.S. Petrich, P.R. Pukite, G.J. Whaley, A.S. Arrott, Surf. Sci. 216 (1989) 222].Method of solution:Epitaxial growth of thin films is modeled by a set of non-linear differential equations [P.I. Cohen, G.S. Petrich, P.R. Pukite, G.J. Whaley, A.S. Arrott, Surf. Sci. 216 (1989) 222]. The Runge-Kutta method with adaptive stepsize control was used for solving initial value problem for non-linear differential equations [W.H. Press, B.P. Flannery, S.A. Teukolsky, W.T. Vetterling, Numerical Recipes in Pascal: The Art of Scientific Computing; first ed., Cambridge University Press, 1989; See also: Numerical Recipes in C++, second ed., Cambridge University Press, 1992].Typical running time: The typical running time is machine and user-parameters dependent.Unusual features of the program: The program is distributed in the form of a main project Growth.dpr file and an independent Rhd.pas file and should be compiled using Object Pascal compilers, including Borland Delphi.     

Constrained MC for QCD evolution with rapidity ordering and minimum kT

With the imminent start of LHC experiments, development of phenomenological tools, and in particular the Monte Carlo programs and algorithms, becomes urgent. A new algorithm for the generation of a parton shower initiated by the single initial hadron beam is presented. The new algorithm is of the class of the so-called “constrained MC” type algorithm (an alternative to the backward evolution MC algorithm), in which the energy and the type of the parton at the end of the parton shower are constrained (predefined). The complete kinematics configurations with explicitly constructed four momenta are generated and tested. Evolution time is identical with rapidity and minimum transverse momentum is used as an infrared cut-off. All terms of the leading-logarithmic approximation in the DGLAP evolution are properly accounted for. In addition, the essential improvements towards the so-called CCFM/BFKL models are also properly implemented. The resulting parton distributions are cross-checked up to the 10⁻³ precision level with the help of a multitude of comparisons with other MC and non-MC programs. We regard these tests as an important asset to be exploited at the time when the presented MC will enter as a building block in a larger MC program for W/Z production process at LHC.     

Photofragment image analysis using the Onion-Peeling Algorithm

With the growing popularity of the velocity map imaging technique, a need for the analysis of photoion and photoelectron images arose. Here, a computer program is presented that allows for the analysis of cylindrically symmetric images. It permits the inversion of the projection of the 3D charged particle distribution using the Onion Peeling Algorithm. Further analysis includes the determination of radial and angular distributions, from which velocity distributions and spatial anisotropy parameters are obtained. Identification and quantification of the different photolysis channels is therefore straightforward. In addition, the program features geometry correction, centering, and multi-Gaussian fitting routines, as well as a user-friendly graphical interface and the possibility of generating synthetic images using either the fitted or user-defined parameters.

Program summary

Title of program: Glass OnionCatalogue identifier: ADRYProgram Summary URL:http://cpc.cs.qub.ac.uk/summaries/ADRYProgram obtainable from: CPC Program Library, Queen's University of Belfast, N. IrelandLicensing provisions: noneComputer: IBM PCOperating system under which the program has been tested: Windows 98, Windows 2000, Windows NTProgramming language used: Delphi 4.0Memory required to execute with typical data: 18 MwordsNo. of bits in a word: 32No. of bytes in distributed program, including test data, etc.: 9 911 434Distribution format: zip fileKeywords: Photofragment image, onion peeling, anisotropy parametersNature of physical problem: Information about velocity and angular distributions of photofragments is the basis on which the analysis of the photolysis process resides. Reconstructing the three-dimensional distribution from the photofragment image is the first step, further processing involving angular and radial integration of the inverted image to obtain velocity and angular distributions. Provisions have to be made to correct for slight distortions of the image, and to verify the accuracy of the analysis process.Method of solution: The “Onion Peeling” algorithm described by Helm [Rev. Sci. Instrum. 67 (6) (1996)] is used to perform the image reconstruction. Angular integration with a subsequent multi-Gaussian fit supplies information about the velocity distribution of the photofragments, whereas radial integration with subsequent expansion of the angular distributions over Legendre Polynomials gives the spatial anisotropy parameters. Fitting algorithms have been developed to centre the image and to correct for image distortion.Restrictions on the complexity of the problem: The maximum image size (1280×1280) and resolution (16 bit) are restricted by available memory and can be changed in the source code. Initial centre coordinates within 5 pixels may be required for the correction and the centering algorithm to converge. Peaks on the velocity profile separated by less then the peak width may not be deconvolved. In the charged particle image reconstruction, it is assumed that the kinetic energy released in the dissociation process is small compared to the energy acquired in the electric field. For the fitting parameters to be physically meaningful, cylindrical symmetry of the image has to be assumed but the actual inversion algorithm is stable to distortions of such symmetry in experimental images.Typical running time: The analysis procedure can be divided into three parts: inversion, fitting, and geometry correction. The inversion time grows approx. as R³, where R is the radius of the region of interest: for R=200 pixels it is less than a minute, for R=400 pixels less then 6 min on a 400 MHz IBM personal computer. The time for the velocity fitting procedure to converge depends strongly on the number of peaks in the velocity profile and the convergence criterion. It ranges between less then a second for simple curves and a few minutes for profiles with up to twenty peaks. The time taken for the image correction scales as R² and depends on the curve profile. It is on the order of a few minutes for images with R=500 pixels.Unusual features of the program: Our centering and image correction algorithm is based on Fourier analysis of the radial distribution to insure the sharpest velocity profile and is insensitive to an uneven intensity distribution. There exists an angular averaging option to stabilize the inversion algorithm and not to loose the resolution at the same time.     

L1PMA: A Fortran 77 package for best L1 piecewise monotonic data smoothing

Fortran 77 software is presented for the calculation of a best L₁ approximation to n measurements that include random errors by requiring k−1 sign changes in the first divided differences of the approximation or equivalently k monotonic sections, alternately increasing and decreasing. A dynamic programming algorithm separates the measurements into optimal disjoint sections of adjacent data and applies to each section a single L₁ monotonic calculation. The most distinctive feature of the algorithm is that it terminates at a global minimum in at most n³+O(kn²) computer operations, although this calculation can exhibit O(n^k) local minima, because the optimal positions of the turning points are also unknowns of the optimization process. The arithmetic operations involved in this calculation are comparisons mainly spent in finding the medians of subranges of data during the monotonic calculations. The package employs techniques for median and for best L₁ monotonic approximation, while full details of these techniques are specified. The package has been applied and tested on a variety of data that have substantial differences and showed quadratic behaviour in n. Some numerical results demonstrate the performance of the method. Further, there is a commentary on the division of the code into subroutines. Driver programs and numerical examples with output are provided to help new users of the method. Besides that piecewise monotonicity is a property of a wide range of functions, an important application of the method is in estimating turning points of a function from some noisy measurements of its values.