The TDHF Code Sky3D

The nuclear mean-field model based on Skyrme forces or related density functionals has found wide-spread application to the description of nuclear ground states, collective vibrational excitations, and heavy-ion collisions. The code Sky3D solves the static or dynamic equations in a three-dimensional Cartesian mesh with isolated or periodic boundary conditions and no further symmetry assumptions. Pairing can be included in the BCS approximation. The code is implemented with a view to allow easy modifications for including additional physics or special analysis of the results.


External routines/libraries: LAPACK, FFTW3
Nature of problem: The time-dependent Hartree-Fock equations can be used to simulate nuclear vibrations and collisions between nuclei for low energies. This code implements the equations based on a Skyrme energy functional and also allows the determination of the ground-state structure of nuclei through the static version of the equations. For the case of vibrations the principal aim is to calculate the excitation spectra by Fourier-analyzing the time dependence of suitable observables. In collisions, the formation of a neck between nuclei, the dissipation of energy from collective motion, processes like charge transfer and the approach to fusion are of principal interest.
Solution method: The nucleonic wave function spinors are represented on a three-dimensional Cartesian mesh with no further symmetry restrictions. The boundary conditions are always periodic for the wave functions, while the Coulomb potential can also be calculated for an isolated charge distribution. All spatial derivatives are evaluated using the finite Fourier transform method. The code solves the static Hartree-Fock equations with a damped gradient iteration method and the time-dependent Hartree-Fock equations with an expansion of the time-development operator. Any number of initial nuclei can be placed into the mesh in with arbitrary positions and initial velocities.
Restrictions: The reliability of the mean-field approximation is lim-

Introduction
The vast majority of microscopic models of many-body systems rely on a description in terms of the single-particle (s.p.) wave functions. Among them, self-consistent mean-field models (SCMF) automatically generate the optimal one-body potentials corresponding to the s.p. wave functions. A rigorous SCMF is the Hartree-Fock theory (HF) where the s.p. wave functions are determined variationally for a given two-body interaction [1,2]. A more practical approach is provided by the Density Functional Theory (DFT), which incorporates the involved many-body effects into effective interactions, or effective energy-density functionals. This is a very efficient and successful scheme, widely used in electronic systems [3]. Straightforward HF is unsuitable for nuclei because the free-space twonucleon force contains a strong short-range repulsion requiring renormalization in the nuclear medium. For this reason, nuclear SCMFs necessarily employ effective interactions or functionals although they often carry the label HF as, e.g, in the Skyrme Hartree-Fock (SHF) method. There are relativistic as well as non-relativistic approaches. For an extensive review see [4].
The description of dynamical processes is even more demanding than the modeling of structure. SCMFs are the also the first method of choice in this domain. The natural extension of HF is time-dependent HF (TDHF) which was proposed as early as 1930 in [5]. Earlier applications were restricted to the linearized regime covering small amplitude motion, see, e. g., [6]. Large scale TDHF calculations became possible in the last few decades with the increasing computing capacities. Again, as in the static case, true TDHF calculations make sense only for electronic systems and even there they are very rare. The overwhelming majority of dynamical SCMF calculations employ, in fact, time-dependent DFT (TDDFT). In electronic systems, this amounts to the time-dependent local density approximation (TDLDA) [3], which is widely used in atoms, molecules, and solids; for examples in nanoparticles see, e. g., [7]. Dynamical SCMFs in nuclei also stay at the level of TDLDA even if they are often named TDHF which happens particularly for dynamical calculations using the Skyrme energy functional. Nuclear TDHF started about forty years ago [8] and has developed since then into a powerful and versatile tool for simulating a great variety of dynamical scenarios. Earlier applications were based mainly on non-relativistic TDHF using the effective Skyrme energy functional [9,10]. Due to higher numerical demands, relativistic calculations appeared somewhat later [11], but have developed meanwhile equally well to a widely used tool [12,13].
In this paper, we present a code for TDHF calculations on the basis of the non-relativistic Skyrme energy functional. The code uses a fully three dimensional (3D) representation of wave functions and fields on a Cartesian grid in coordinate space. There are no symmetry restrictions and the full Skyrme energy functional is used including the spin-orbit and most important time-odd terms. Such fully-fledged 3D calculations became possible only over the last decade with the steadily increasing computing capabilities. In fact, early TDHF studies all used restricted representations, axial symmetry and/or reflection symmetries. This limited the possible applications. TDHF experienced a revival during the last ten years when unrestricted 3D calculations became possible. There are several groups performing large scale TDHF studies for various scenarios of nuclear dynamics, see, e.g., [14,15,16,17,18,19,20]. For a recent review see [21]. Such calculations have clearly outgrown the developmental stage. It is an appropriate time to give a broader public access to a 3D TDHF code. This is the goal of this paper. Skyrme HF covers such a broad range of physical phenomena and is relatively involved that efficient computational treatment of 3D simulations requires elaborate numerical methods. We shall make an effort to explain the many necessary ingredients in a comprehensive, and yet compact, manner.
Most recent Skyrme density functionals contain terms such as fractional powers of the density that cannot be related to a two-or three-body interaction. In that sense, the present code solves the TDDFT rather than TDHF equations. Nevertheless, we prefer to keep the name TDHF since it is associated historically with this large field of nuclear reaction theory.

Intended applications
The code Sky3D solves the static Hartree-Fock as well as the time-dependent Hartree-Fock (TDHF) equations for interactions of Skyrme-force type in a general three-dimensional geometry. No symmetries of any kind are assumed, so that the code can be used for a wide variety of applications in nuclear structure, collective excitations, and nuclear reactions; of course within the limitations of mean-field theory.

Specific model implemented
The code in the presented version contains a useful selection of terms in the Skyrme force but by no means all terms that have been included in some recent works. It should still be useful, because (1) for many interesting applications the interest is semi-quantitative so that a Skyrme force fitted with the latest models is not necessary -usually a selection of forces is desired to look at the variability of results, but not a high-accuracy fit of data; (2) the code is written in such a way that additional terms can be added easily. The coding corresponds one-by-one to the analytic formulas in most places except where efficiency demands reordering the calculations.

The single-particle basis
In a mean field theory one seeks to describe the many-body system exclusively in terms of a set of single-particle wave functions ψ α with fractional occupation amplitudes v α , i.e.
{ψ α , v α , α = 1, ..., Ω} where Ω denotes the size of the active s.p. space. The occupation amplitude can take values continuously in the interval [0, 1]. The complementary non-occupation amplitude is u α = 1 − v 2 α . A formal definition of the BCS mean-field state reads |Φ = α>0 u α + v αâ where |0 is the vacuum state,â + α the generator of a Fermion in state ψ α , andᾱ the time reverse partner to state α. We will use variation of the BCS amplitudes v α only in the static part of even-even nuclei where the time reverse partner is unambiguously defined. These amplitudes are kept frozen during dynamical propagation whenever BCS pairing is used for a nucleus.

Local densities and currents
The Skyrme-energy-density functional is defined in terms of only a few local densities and currents. These are the timeeven fields the time-odd fields plus the exchange term in the Slater approximation [23]. The formula is • E pair : the pairing energy. It consists of a contact pairing interaction involving the pairing densities ξ q augmented by an optional density dependence. The formula is It contains a continuous switch, the parameter ρ 0,pair . A pure δ-interaction (DI), also called volume pairing, is recovered for ρ 0,pair −→ ∞. The general case is the density dependent δ-interaction (DDDI). A typical value near matter equilibrium density ρ 0,pair = 0.16 fm −3 concentrates pairing to the surface. The most flexible choice is to consider ρ 0,pair as an additional free parameter. Actual adjustments with this option deliver a form of the pairing functional which stays in between the extremes of volume and surface pairing [24].
The term E corr stands for all additional corrections from correlations beyond the mean field that might be added. Most calculations include at least the center-of-mass correction E cm . For deformed nuclei this should be augmented by a rotational correction and for soft nuclei by correlations from all low-energy quadrupole motions [25]. So far, these correlations are usually added a posteriori after static calculations. A fully variational treatment and a dynamical propagation of the c.m. correction is extremely involved and usually not considered. We aim predominantly at dynamical simulations and thus here ignore all such correlations, using This has to be taken into account when comparing ground state energies from this code with results from other static codes which contain a c.m. correction E cm 0. However, one has to take into account that there are two different ways to perform the c.m. correction. They are selected by the parameter zpe. One strategy is to augment the energy by E cm a posteriori. We associate this option with zpe=1. This means that in our code nothing is done. The other strategy is to modify the nucleon mass by m −→ m − m/A and to not have a posteriori correction. This way is chosen in a couple of traditional parametrizations, e.g., in SkM * [26]. We keep this option for consistency and associate it with zpe=1. But note that this choice runs into inconsistencies in collisions and fragmentation as it employs only the total mass number A and cannot account for the masses of the fragments.
The functional in the above form contains the minimal number of terms which are needed to guarantee Galilean invariance [27,22] and so to allow performance of TDHF calculations which respect all basic conservation laws. We ignore the tensor spin-orbit terms and spin-spin couplings [4,28,22]. These may be important for magnetic excitations [22] and odd nuclei [29] which are, however, not the prime focus of TDHF studies.

Force coefficients
The above formulation in terms of the Skyrme energy functional introduces the force parameters b 0 , b ′ 0 , ... b ′ 4 naturally as the factors in front of each contribution in the terms (5c). Traditionally, the functional is deduced from a Skyrme force which is a density-dependent, zero-range interaction [30]. The t and x coefficients in this Skyrme-force definition are related to the b coefficients in the functional definition as The coefficient b ′ 4 is usually fixed as b ′ 4 = 1 2 t 4 for most traditional Skyrme forces. More recent variants of Skyrme forces (SkI3 etc.) handle it as a separate free parameter [31]. In addition to the b and b ′ parameters, there is the power coefficient in the (originally) three-body term, which is usually called α, but in the code is referred to as power. The input of the force to the code is done in terms of the t i , x i parameters.

The single-particle Hamiltonian
The mean-field Hamiltonianĥ is derived from the energy functional of Section 2.2.3 by variation ∂ ψ * α E =ĥψ α . It reads in detailĥ with q ∈ {p, n} as usual distinguishing the different Hamiltonians for protons and neutrons. For the protons the Coulomb potential is also added. The first term is the local part of the mean field, which acts on the wave functions like a local potential. It is defined as Next comes the "effective mass", which replaces the standard 2 2m factor by the isospin and space-dependent expression Note that the Skyrme force definitions contain the first term (nucleon mass) as a parameter which varies slightly from parametrization to parametrization and may be different for protons and neutrons. The spin-orbit potential is The above three terms involve the time-even densities. Dynamical effects come into play with the next terms which include the time-odd contributions from current and spin-density:

Coupling to external fields
For the dynamic case, the system can also be coupled to an external excitation field, to study collective response such as in giant resonances. The present code only implements a very simple case, since it is expected that most serious applications will need modifications, which are quite easy to incorporate.
The external field is introduced as a time-dependent, local operator where f (t) carries the temporal profile of the excitation mechanism, F q ( r) is some local operator, and η tunes the overall strength. The spatial distribution F q ( r), is allowed to be different for the two isospins q. Typical examples are isoscalar and isovector multipole operators as, e.g., the isoscalar quadrupole The prefactor η is a strength parameter which allows scanning different excitation strengths easily while keeping the temporal and spatial profiles the same. It should be noted that the absolute magnitude of the perturbing potential by itself usually has little direct meaning. What counts is the excitation energy caused by the perturbation and subsequently the magnitude of vibrations in the observables (such as the time-dependent quadrupole moment). An exception is, e. g., the simulation of the close approach of another nucleus that stays external to the computational grid, where the potential is uniquely defined.
One important point remains to be noted concerning the spatial profile F q ( r). This can be illustrated by the quadrupole operator. Let us assume an instant where A > 0 and f (t) > 0. For then the operator ∝ 2z 2 − x 2 − y 2 leads to a perturbing potential U ext which is binding in z-direction but asymptotically unbound in the x-and y-directions. This can cause unphysical effects in case of large strengths and/or numerical boxes. For this reason it is useful to have a cut-off by a Woods-Saxon like function according to [32] where r 0 and ∆r are parameters describing a transition region sufficiently outside the nucleus, but also sufficiently small to maintain binding. Another problem associated with the external field is that in general it will not be periodic but instead have discontinuities on the boundary when crossing into the neighboring cell. If damping is sufficiently strong, the field may be practically zero on the boundary and thus become periodic. Another solution for this problem is to make the field explicitly periodic by replacing the coordinates with periodic substitutes. The exact formulation depends on the specific field used; for the abovementioned quadrupole operator, which depends only on the squares of the coordinates, e. g., substituting with x L = nx ∆x the period interval, will provide the proper behavior as the sine squared has a period of π and there is no unphysical decrease of this function near the boundary. Of course the analogous transformation has to be applied to y and z.
The time dependence f (t) is modeled as a short pulse centered around some excitation frequency ω. The code allows a choice between two pulse envelopes:

A Gaussian of the form
with peak time τ 0 and width ∆τ.
Broad envelopes provide a high frequency resolution and so concentrate the excitation around the driving frequency ω. Short pulses lose frequency selectivity and excite a broad band of frequencies.
The extreme case is an infinitely short pulse, ∆τ −→ 0. It amounts eventually to an instantaneous boost of the initial wave functions which can be expressed as a phase factor according to where ψ k,0 stands for the stationary wave function before boost. This instantaneous boost, being infinitely short, is insensitive to the problem of asymptotically unstable potentials and allows the use of a driving field F q without damping (9b) which, in turn, simplifies spectral analysis (see Section 2.5.3).
The effect of the boost (10) can be understood by virtue of the Madelung decomposition of a complex wave function φ( r) = χ( r) exp(iS ( r)) with χ and S being real. A straightforward calculation leads to the probability flow density as  = χ 2 ∇S /m. Dividing by the density χ 2 then produces the "probabilityflow velocity" v = ∇S /m. An illustrative example is the plane wave where χ =constant and S = k · r which yields the correct result v = k/m. In the present case, assuming that the static wave functions themselves can be written as real functions, we get an initial velocity v = −∇F q , i. e., just in the direction of the classical force resulting from the "velocity potential" F q (r).

Static Hartree-Fock 2.4.1. The coupled mean-field and BCS equations
The stationary equations are obtained variationally. Variation with respect to the single-particle wave functions ψ α yields the mean field equations [33,2] hψ α = ε α ψ α , whereĥ is the mean-field Hamiltonian (8a) and ε α is the singleparticle energy of state α. Simultaneous variation of ψ α together with the occupation amplitude v α yields the Hartree-Fock-Bogolyubov equations [33,34,4] which complicate Eq. (11a) by non-diagonal terms on the right-hand-side [35]. We employ here the BCS approximation which exploits time-reversal symmetry where each single-particle state has a time reversed partner ψ α ↔ ψ α as was already implied in the pairing densities ξ q in Eq. (4). Each pair of time-conjugate states is associated with an occupation amplitude v α . These are determined by the BCS equation ε α − ǫ F,q α u 2 α − v 2 α = ∆ α u α v α which can be resolved to a closed expression for the occupation amplitudes as q α denotes the nucleon type to which state α belongs, α ∈ q means all states of type q, and N q is the number of nucleons of type q (identified as N p = Z and N n = N). The Fermi energies ǫ F,q α serve to regulate the average particle number to the required values N q . Here, the space of pairing-active states is just the space of states actually included in the calculation. The results of BCS pairing depend slightly on the size of the active space [33,34]. We recommend using about q single-nucleon states, which comes closest to the dynamical pairing space of Ref. [36].

Iterative solution
The coupled mean-field and BCS equations (11) are solved iteratively. The wave functions are iterated with a gradient step which is accelerated by kinetic-energy damping [37,38] whereT =p 2 /(2m) is the operator of kinetic energy, O means orthonormalization of the whole set of new wave functions, and the upper index indicates the iteration number. Note that this sort of kinetic-energy damping is particularly suited for the fast Fourier techniques that we use in the present code. The damped gradient step has two numerical parameters, the step size δ and the damping regulator E 0 . The latter should be chosen typically of the order of the depth of the local potential U q . In practice, we find E 0 = 100 MeV a safe choice. The step size is of order of δ = 0.1...0.8. Larger values yield faster iteration, but can run more easily into pathological conditions. The optimum choice depends somewhat on the Skyrme parametrization. Those with effective mass m * /m ≈ 1 allow larger δ values. Low m * /m may require a reduction in the step size. After each such wave function step, the BCS equations (11b-11d) are solved with ε α = ψ α |ĥ|ψ α , the densities are updated (Eqs. [2][3][4], and new mean fields computed. (8). This then provides the starting point for the next iteration. The process is continued until sufficient convergence is achieved. We consider as the convergence criterion the average energy variance, or fluctuation, of the single particle states (13c) The single particle energy ε α is defined here as an expectation value. It finally becomes an eigenvalue in Eq. (11a) if the iteration has converged to ∆ε α ≈ 0. Vanishing total variance ∆ε signals that we have reached minimum energy, i.e. a solution of the mean-field plus BCS equations. One has to be aware, however, that this may be only a local minimum (isomeric state). It requires experience to judge whether one has found the absolute energy minimum. In case of doubt, one should redo a couple of static iterations from very different initial configurations. This raises the question of how to initialize the iteration. We take as a starting point the wave functions of the deformed harmonic oscillator (see point 1 in Section 2.8). These are characterized by n = (n x , n y , n z ), the number of nodes in each direction. We stack the wave functions in order of increasing oscillator energy ǫ (0) α = ω x n x + ω y n y + ω z n z and stop if the desired number of states is reached. The deformation of the initializing oscillator influences the initial state in two ways: first, through the deformation of the basis wave functions as such, and second, through the energy ordering of the ǫ (0) α and corresponding sequence of levels built. Variation of initial conditions means basically a variation of the oscillator deformation. For example, the iteration will most probably end up in a prolate minimum if the initial state was sufficiently prolate, and in an oblate minimum after an oblate initial state. It depends on the nucleus which one is the absolute minimum.

The time-dependent mean-field equations
The TDHF equations are determined from the variation of the action with respect to the wave functions ψ + α where the energy is given as in Eq. (5) [2]. This yields the TDHF equation whereĥ is, again, the mean-field Hamiltonian (8a). For simplicity, we are not considering variation of the occupation amplitude in the time-dependent case. The occupation amplitudes obtained from static iteration are kept frozen during time evolution. For studies of mean-field flow at moderate excitations (heavy-ion collisions, giant resonances) this approximation is legitimate. However, a study of truly low energy dynamics in the range of a few MeV (soft vibrations, fission) requires a full time-dependent Hartree-Fock-Bogolyubov treatment and should not be undertaken with the present code.

Time development algorithm
The TDHF equation (14) can be formally resolved into an integral equation as whereÛ is the time-evolution operator andT the time-ordering operation. This time evolution is unitary, thus conserving orthonormalization of the single-particle wave functions, and it conserves the total energy (5) provided that there are no timedependent external fields. To convert this involved operator into an efficiently computable but also sufficiently accurate form a predictor-corrector strategy is used: 1. In a first step (predictor), we determine the single-particle Hamiltonian at midtimeĥ(t+∆t/2). To that end, a trial step by ∆tψ is performed using the mean fieldĥ(t) at initial time t. This produces an estimateh of the mean field Hamiltonian at t + δt. The Hamiltonian ad midtime is then obtained as the averageĥ(t + δt) = ĥ (t) +h /2. 2. In a second step (corrector), the trial-step wave functions {ψ α } are used to compute an estimate (predictor) of the mean-field at the final timeĥ pre in the standard manner by accumulating the densities and composing these to the Skyrme mean-field Hamiltonian (8). The predictor for the Hamiltonian at midtime t + ∆t/2 then becomesh ≈ ĥ (t) +ĥ pre /2. This is used to finally perform a full time step (again with frozen Hamiltonian, but nowh) 3. In both cases the operator exponential is evaluated by a Taylor series expansion up to order m: whereĥ is the actual mean field in step (16), or (17) respectively.ĥ n ψ is computed in straightforward manner by successive application of the mean field Hamiltonian, i.e.ĥ n ψ =ĥ(...(ĥ ntimes ψ)...).
The Taylor expansion spoils strict unitarity of the exponential exp − i ĥ ∆t and energy conservation. We turn this flaw into an advantage and use norm conservation as well as energy conservation (if it applies) as counter-check of the quality of the step along the propagation. The reliability depends, of course, on a proper choice of the numerical parameters in this step which are the step size ∆t and the order of the Taylor expansion m. The step size is limited by the maximum possible kinetic energy and by the typical time scales of changes in the mean field h. The maximum kinetic energy, in turn, depends on the grid spacing as ∝ ∆x −2 . A choice of ∆t = 0.1 − 0.2fm/c is applicable in connection with ∆x = 0.7 − 1/fm. For the order of Taylor expansion, one needs at least m = 4. One may also consider m = 6. Higher orders m are not necessary for the typical values of ∆t.

Collective excitations
Giant resonances are prominent excitation modes of nuclei. Best known is probably the isovector giant dipole resonance, but there are many more modes depending on isospin and angular momentum. The typical resonance energies lie in a region from 10 to 30 MeV where the present TDHF code with frozen occupation numbers is applicable because the relevant energy range lies far above the pairing gap (1-2 MeV). The generation of these modes is particularly simple within the present TDHF treatment. One first produces a stationary state as outlined in Section 2.4 and then triggers the excitation by a timedependent external field as described in Section 2.3. A broad pulse allows triggering particular excitation energies. An infinitely short pulse amounts to an instantaneous boost.
The boost is a generic excitation of a system which gives the same weight to all frequencies. It is thus ideally suited for analyzing in an unbiased manner the excitation spectra of a system. This, in turn, allows a thorough spectral analysis. To obtain the spectral distribution of isovector dipole strength, one applies a boost with small strength η and F q =D ∝ r 1 Y 10 τ z the isovector dipole operator. The Slater determinant |Φ(t) is propagated in TDHF for a sufficiently long time while recording the dipole moment D(t) = Φ(t)|D|Φ(t) . The dipole strength is finally extracted from the Fourier transformD(ω) as S D (ω) = ℑ D (ω) /η. The straightforward Fourier transform leads to artifacts if the dipole signal has not fully died out at the end of the simulation time. In the general case, some filtering is necessary to suppress artifacts from cutting the signal at a certain final time [39]. In practice, it is most convenient to use filtering in the time domain by damping the signal D(t) towards the final time. A robust choice is (19) where t final is the final time of the simulation. This guarantees that the effective signal D fil vanishes at the end of the interval. The cos n profile switches off gently and leaves as much as possible from the relevant signal at early times. The parameter n fil determines the strengths of filtering. Value of order of 4-6 are recommended to suppress the artifacts safely. For a detailed description of this spectral analysis see [40]. For typical applications in nuclear physics see [41]. It is to be noted that the code does not include this final step of spectral analysis. The time dependent signals are printed on the protocol files monopoles.res, dipoles.res, and quadrupoles.res. It is left to the user to perform the final steps towards a spectral distribution. A word is in order about t final . It determines the resolution of the spectral analysis. The corresponding energy bins are given by δE exc = π/t final . Windowing effectively reduces the time span in which relevant information is contained and roughly doubles the relevant δE exc . For example, to obtain a spectral resolution of 1 MeV, one needs to simulate up to about 1200 fm/c.
Although excitation spectra are one of the most basic properties of the system, there are many other dynamical features of interest. The multipole signals in the time domain (printed in the protocol files) are as such interesting quantities. One can have, e.g., a look at cross-talk between the multipole channels. It is particularly interesting to study excitation dynamics for varying excitation strength η, from the regime of linear response (small η) deep into the non-linear regime. It is inefficient to perform a full three-dimensional TDHF calculation to obtain linear-regime excitation spectra for spherical nuclei. This is better done in a dedicated RPA calculation on a spherical basis (see, e.g., [42]) for which there exist an overwhelming multitude of codes. The realm of TDHF calculations of nuclear excitations are spectra in deformed systems, stability analysis of exotic configurations, and in particular non-linear dynamics.
There are many more details worth looking at. One may check the densities and currents to visualize the flow pattern associated with a mode. A most elaborate analysis deals with a phase-space picture of nuclear dynamics by virtue of the Wigner transformation [43]. The code allows saving all ingredients needed for such elaborate analysis in dedicated output files, see Section 7.5. It is left to the user to work out the further steps of the analysis.

Nuclear reactions
Collisions of nuclei are a prime application of nuclear TDHF. They were, in fact, the major motivation for its realization [44,10]. The present code is designed to initialize such collision scenarios in a most flexible manner. We start by explaining the simplest case of a collision of two nuclei. First, we prepare the ground states of the two nuclei as explained in Section 2.4. where I = 1, 2 labels the two nuclei are shifted to new centers R I where the distance | R 2 − R 1 | should be sufficiently large that the nuclei have negligible overlap and negligible Coulomb distortion from the other nucleus (the latter condition usually only loosely fulfilled). The shifted wave functions ψ α,I ( r − R I , s) (stat) are obtained by interpolation on the grid. It is obvious that the collisions need a larger numerical box than the static Hartree-Fock calculations. Thus, we may compute the static wave functions on a smaller box since we are shifting and interpolating the result anyway for dynamical initialization. At this point we have the nuclei resting at a safe distance. To set them in motion we need to give each nucleus a momentum P 1 = − P 2 (note that the total momentum of the combined system still vanishes). Consequently, the initial configuration is given by the Slater state built from the shifted and boosted single-particle wave functions (see fragment initialization, point 2 in Section 2.8) The distance between the nuclei is large, but inevitably finite. This may induce minor violations of orthonormality. Thus the full set of wave functions (20) is orthonormalized as a final step of initial preparation.
The occupation amplitudes v α,I are taken over from the static solution and frozen along the dynamical evolution. In fact, most of the collision studies will be principally to explore the dynamical features in the regimes of fusion and inelastic collisions. It is then recommended to use the conceptually simplest and most robust strategy, namely to start from simple Slater states (not BCS states) for the two nuclei. This means we that in most cases the static solution is calculated without pairing, fixing v α = 1 and including just as many states as there are nucleons.
The above example deals with two initial fragments. The code is more flexible than that. It allows an initial state composed from several fragments. The strategy remains the same as for the binary system. It is simply repeated for each new fragment.
The time evolution is performed as outlined in Section 2.5.2. It requires some effort to visualize the complex dynamics which emerges in collisions. A rough picture is, again, provided by the multipole moments. The quadrupole moment, e.g., can serve as a measure of stretching of the total system. Small values indicate a compound nucleus while asymptotically growing values signal fragmentation. One may want, particularly in case of collisions, more detailed pictures of the flow as, e.g., snapshots of the density of current distributions and, ultimately, a full phase space picture [43]. Material for that can be output on demand, as detailed in Section 7.5. Again, we leave it to the user to extract the wanted information and to prepare it for visualization.

Observables
It was already mentioned in Sections 2.5.3 and 2.5.4 which observables may be used to analyze nuclear dynamics. We here briefly summarize the observables computed and output in the code and indicate how further observables may be extracted. Basic features of the description by the Skyrme energy-density functional are, of course, energy and densities.

Multipole moments
The gross features of the density distribution are well characterized by the multipole moments. The most important moment is the center of mass (c.m.) where A = d 3 r ρ( r) is the total mass number and "type" can refer to proton c.m. from ρ p , neutron c.m. from ρ n , isoscalar or total c.m. from the total density ρ = ρ p +ρ n ≡ ρ T =0 , or isovector moment related to the isovector density The definition of R (type) directly employs the Cartesian coordinate r i . The same holds for the Cartesian quadrupole tensor again for the various types as discussed above. The matrix Q kl is not invariant under rotations of the coordinate frame. There is a preferred coordinate system: the system of principle axes. It is obtained by diagonalizing Q kl . The quadrupole matrix in the principle-axis frame thus has only three non-vanishing entries Q xx , Q yy , and Q zz together with the trace condition Q xx + Q yy +Q zz = 0. Although straightforward to define, higher Cartesian moments are messy and inconvenient to use. The spherical multipole moments are deduced from where r = | r| and Y lm are spherical harmonics. The most important ones are monopole (l = 0), dipole (l = 1), and quadrupole (l = 2) moments. The latter are often expressed as dimensionless quadrupole moments with R = r 0 A 1/3 a fixed radius derived from the total mass number A. This, again, could be defined for any "type", but is used, in practice, mainly for the isoscalar moments. r rms is the rootmean-square radius of the total density. The r.m.s. radii are defined as where "type" can be proton, neutron, or total. The isovector variant does not make sense here. The dimensionless moments have the advantage of being free of an overall scale which was removed by the denominator AR 2 . They allow characterization of the shape of the nucleus. However, the general a m are not invariant under rotations of the coordinate frame. We obtain a unique characterization by transforming to the principle-axis system. These are defined by the conditions a ±1 = 0 and a 2 = a −2 . There remain only two shape parameters a 0 and a 2 . These are often reexpressed as total deformation β and triaxiality γ, often called Bohr-Mottelson parameters, through Triaxiality γ is handled like an angle. It can, in principle, take all values between 0 o and 360 o , but physically relevant parameters stay in the 0 . . . 60 • range. The other sectors correspond to equivalent configurations [33]. We supply printouts of all the above variants of multipole moments to allow a most flexible analysis.

Alternative way to evaluate the total energy
The key observable is the total energy E tot . It is computed as given in Eqs. (5). More detailed energy observables are provided by the s.p. energies (13c). These can also be used to compute the total energy. The traditional HF scheme deals with pure two-body interactions and exploits that to simplify [33] where t α = ψ α |T |ψ α is the s.p. kinetic energy. This is possible because where u α is the s.p. mean-field potential energy and v the twobody interaction, and The Skyrme force does not simply have this two-body structure. Still the total energy is very often computed along the strategy of Eq. (22). However, the density dependence requires augmenting this recipe by a rearrangement energy which accounts for a contribution missing in the simple recipe (22). The extension to the Skyrme energy thus reads In the code the total energy is computed both ways, from the straightforward Skyrme energy (5) as well as from the above recipe (23). Numerically these values are close but not identical.

Data types
In the module Params a type db is defined for 12-digit accuracy and on all present machines should amount to double precision, i. e., REAL(db) and COMPLEX(db) are actually REAL (8) and COMPLEX(8) for most compilers. Keeping the symbolic type throughout of course makes the code more flexible for future hardware. Using single precision is not recommended.
Since the external libraries are based on C or Fortran-77 coding, special care has to be taken in this respect. The FFTW3 library stores its plans in 8-byte integers, which in the modules Coulomb and Fourier are defined using the type C LONG from the system-supplied module ISO C BINDING. If this is not available, they may be defined as INTEGER (8) or if that also causes problems, DOUBLE PRECISION, which certainly corresponds to at least 8 bytes. For the LAPACK routines, doubleprecision real and complex variables are necessary, which we also define using "db''. If db should ever be changed in such a way that REAL(db) no longer corresponds to 8 bytes, different LAPACK routines must be selected.
Note that data conversion needs some care. If a and b are double precision, CMPLX(a,b) is not; according to the standard it returns the default accuracy, which on many machines will still be single precision. Thus that for example in EXP(CMPLX(a,b)) the exponential would be evaluated in single precision. That is why in the code the expression CMPLX(a,b,db) is used consistently; this is also safe for future changes in accuracy.
The other conversion functions used are safe. AIMAG is generic and REAL reproduces the accuracy of (only) a complex argument.

Grid definition
All wave functions and fields are defined on a three-dimensional regular Cartesian grid of nx by ny by nz grid points. nx, ny, and nz must be even numbers. The physical spacing between the points is given as dx, dy, and dz (in fm). In principle these could be different for the three directions, but since this will lead to a loss of accuracy it is highly recommended to give the same value to all of them. A typical range is 0.5-1.0 fm.
The grid is automatically arranged in such a way that in each direction the same number of grid points are located on both sides of the origin. This means that the three-dimensional origin is in the center of a cubic cell and has the advantage that exact parity properties for the wave functions can be maintained. The coordinate values for e. g., the x-direction are thus: The corresponding values are available in the arrays x(nx), y(ny), and z(nz).

Derivatives
The computation of the kinetic densities and currents and the application of the mean-field Hamiltonian require first and second derivatives at several places in the code. We define them in Fourier space. For simplicity, the strategy is explained here for one dimension. The generalization to 3D is obvious.
The nx discrete grid points x ν in coordinate space are related to the same number of grid points k n in Fourier space (physically equivalent to momentum space) as Note the particular indexing for the k-values. In principle, the values k n = (n − 1)dk for all n are equivalent for the Fourier transform, but for the second half of this range the negative kvalues should be chosen because of their smaller magnitude. For the Fourier expansion, k = −dk and k = (nx − 1)dk are equivalent because of periodicity in k-space.
This complex Fourier representation implies that the function f is periodic with f (x + dx · nx) = f (x). The appropriate integration scheme is the trapezoidal rule which complies with the above summations adding up all terms with equal weight. The derivatives of the exponential basis functions are Computation of the mth derivative thus becomes a trivial multiplication by (ik n ) m in Fourier space. Time critical derivatives are best evaluated in Fourier space using the fast Fourier transformation (FFT). To that end a forward transform (25d) is performed, then the valuesf (k n ) are multiplied by (ik n ) m as given in Eq. (25b) and finally transformed (ik n ) mf (k n ) back to coordinate space by the transformation (25e). This strategy is coded in the subroutines cdervx, cdervy, and cdervz contained in module Levels. It is used for derivatives of wave functions provided the switch TFFT is set. For coding purposes, it is often useful to perform derivatives as a matrix operation directly in coordinate space. The derivative matrices are built by evaluating the double summation of forward and backward transform ahead of time. For the mth derivative this reads From here, the detailed handling depends on the order of derivative. The k n run over the values k n = 0, ±dk, ±2dk, ..., (nx/2 − 1)dk, +dk nx/2. Here the index ordering given in Eq. (25b) does not matter as the index is summed over. Note that the first and the last value come alone while all others come in pairs of ± partners. These pairwise terms can be combined into a sine function for n = 1 and a cosine for n = 2. The derivative matrices thus read in detail A word is in order about the first derivative. The upper point in the k-grid, dk nx/2, is ambiguous. Exploiting periodicity, it could be equally well −dk nx/2. In order, to deal with a ±k symmetric derivative we have anti-symmetrized this last point. The price for this is a slight violation of hermiticity which, however, should be very small as we anyway should not have significant wave-function contributions at the upper edge of the k-grid. The derivative matrices D (m) νν ′ can be prepared ahead of time and are then at disposal for any derivative in the course of the program. Actually, the matrices for the first derivative are prepared in routine sder and for the second derivative in routine sder2, both contained in module Grids. These routines are applied to generate the derivative matrices derv1x, derv1y, derv1z, derv2x, derv2y, and derv2z, for first and second derivatives in the x, y and z directions.
The matrix formulation of the derivatives is used in the code in two ways: on the one hand, the derivatives of the real-valued densities, currents, and mean-field components are always calculated using these derivative matrices, because they are real and the Fourier transform method would require converting them to complex values (using special Fourier techniques for real arrays is in principle possible but has not been worked out yet). In addition the user can switch to using the matrix method everywhere, which may give a slight speed advantage for small grid dimensions.

Boundary conditions
The code uses a periodic Fourier transform to calculates derivatives. This is valid only with periodic boundary conditions. Thus in principle the wave functions and potentials are assumed to be repeated periodically in each Cartesian direction. Because of the short range of the nuclear force, this is not a serious problem in most cases; at higher energies, however, the emission of low-density material from the nuclei can interfere with the dynamics in the neighboring box and cause problems in the conservation of energy and angular momentum; for a detailed discussion see [45]. This is aggravated by the fact that even with periodic boundary conditions periodicity is truly fulfilled only for the wave functions and mean-field components. Since the vector r itself is not periodic but jumps at the boundary, operators such as the orbital angular momentum are not periodic.
On the other hand, for the Coulomb field with its long range this would be clearly wrong. Therefore a computation of the Coulomb potential for the boundary condition of an isolated charge distribution is implemented in addition to the periodic one, see the manual for the module Coulomb. This is selected by the logical input variable periodic and applies only to the Coulomb potential.

Wave function storage
Module Levels handles the single-particle wave functions and associated quantities.
The principal array for the wave function is called psi, which is of type COMPLEX(db). Its dimension is (nx,ny,nz,2,nstloc), where the first three indices naturally refer to the spatial position. The 4th index corresponds to spin: index 1 refers to spin up and 2 to spin down, quantization being along the z-direction.
The last index numbers the wave functions. If the code is run on a single node, the value is nstmax, the total number of single-particle wave functions. They are divided up into neutron and proton states, with the index range given by npmin and npsi. The sub-ranges are: • npmin(1). . .npsi(1) : the neutron states, • npmin(2). . .npsi(2) : the proton states.
If the code is run in parallel (MPI) on several nodes, only nstloc single-particle wave functions are stored on a given node, where nstloc may vary. Pointers are then defined to indicate the relationship between the local index and that in the global array of wave functions. For details see the section on parallelization. There are a number of arrays containing the physical properties of the wave functions, such as the single-particle energy. The names start with sp and they are defined in module Levels. They are not split up in the parallel case, but on each node only the pertinent index positions are used.

Densities and currents
The various densities necessary for constructing the mean field are actually kept in separate arrays and can be output onto data files for later analysis (see subroutine write densities). The dimensioning is (nx,ny,nz,2) with the last index referring to isospin for scalar densities, so that rho(:,:,:,1) is the neutron density and rho(:,:,:,2) the proton density. For vector densities there is an additional index with values 1 to 3 for the Cartesian direction, thus sdens(nx,ny,nz,3,2) containing the spin density in each direction for neutrons and protons.
Since it is often not necessary to keep the neutron and proton contributions separate, subroutine write densities has the option of adding them up before output.

Initialization
A particular strength of the code is its flexible initialization. There are essentially three types of initialization, which can be selected through the input variable nof: 1. Harmonic oscillator: nof=0: this is applicable only to static calculations. The initial wave functions are generated from harmonic oscillator states with initial radii radinx, radiny, and radinz in the three directions. It is advisable to choose the three radii different to avoid being kept in a symmetric configuration for non-spherical nuclei. Note that this is a very simple initialization and has some defects; for example, the initial deformation is controlled more by the occupation of the oscillator states than by the radius parameters. This should eventually be replaced by, e. g., Nilsson wave functions. For this case the type of nucleus is determined by the input numbers nneut and nprot giving the number of neutrons and protons, while npsi can be used to add some unoccupied states (this sometimes leads to faster convergence). 2. Fragment initialization: nof>0: wave functions for a number nof of fragments are read in and positioned in the grid at certain positions. The wave functions are read from files produced by the static code with the file names given by the input filename, they are positioned at centerof-mass positions fcent and given am initial velocity controlled by fboost. The code determines the number of wave functions needed from these data files and also checks the agreement of Skyrme force and grid used. This initialization is used, e.g., for nuclear reactions (see Section 2.5.4).
The number of fragments read in is arbitrary, but there are two special cases: • for nof=1 a single fragment is read in. This can be useful for initializing with static wave functions to study collective vibrations in a nucleus using the TDHF mode. • for nof=2 a special initialization can be done where the initial velocities are not given directly but computed from a center-of-mass energy ecm and an impact parameter b.

User initialization: a user-supplied routine user init
can be employed to set up the wave functions in any desired way. The only condition is that the index ranges etc. are set up correctly and the wave function array psi is filled with the proper values. It was found useful, e. g., to use initial Gaussians distributed in various geometric patterns for α-cluster studies.

Restarting a calculation
Sometimes it is necessary to continue a calculation that was not run to the desired completion because of a machine failure or because the number of iterations or time steps was set too low. In such cases the last wave function file with name wffile, which is generated at regular intervals of mrest iterations or time steps, can be used to initialize a continuation. The program handles this in a simple fashion: if the logical variable trestart is input as TRUE, it sets up an initialization with one fragment (read from the initialization file) placed at the origin and with zero velocity. The only other modifications to the regular setup are then to take the initial iteration number and time from that file instead of starting at zero, as well as suppressing some unneeded initialization steps. This flexible restart makes it possible to use a different grid for the continuation in the sense that the grid spacings must agree, but the new grid can be larger than the old one.

Accuracy considerations
The grid representation and solution methods introduced above depend on several numerical parameters. Their proper choice is crucial for the accuracy and speed of the calculations. In this Section, we want to briefly address the dependence on numerical parameters. An extensive discussion of grid representations and static iteration is found in [38].   Figure 2 shows the sensitivity with respect to the grid spacings ∆x, ∆y, ∆z. The trend is the same for both observables, energy and radius: The results have very high quality and change very little up to ∆x = 0.75 fm. They quickly degrade above that spacing. But even at ∆x = 1 fm, we still find an acceptable quality which suffices for most applications, particularly for large scale explorations. If high accuracy matters, ∆x ≈ 0.75 fm should be chosen; not much is gained by going to even finer gridding. This holds for ground states and moderate excitations. High excitations and fast collisions may require a finer mesh. Note that the maximum representable kinetic energy is E kin,max = ( 2 /2m)(π/∆x) 2 , which amounts to about 200 MeV for ∆x = 1 fm. The actual energies of interest should stay far below this limit. It is an instructive exercise to study uniform center-of-mass motion at various velocities to explore the limits of a given representation.
The number of grid points N x = N y = N z in the tests of Figure 2 were chosen such that the box size was the same in all cases. The actual choice of N x depends sensitively on the system, its size and separation energy. As a rule of thumb, the density decreases asymptotically as ρ ∝ exp (−2 √ 2mε N r/ ) where ε N is the single particle energy of the least bound state. One should aim for at least ρ < 10 −8 at the boundaries. Figure 3 explores the effect of grid spacing for dynamics. Two different spacings are compared for a quadrupole oscillation following an instantaneous quadrupole boost. Practically no difference can be seen for the "safe choice" ∆x = 0.75 fm and the robust choice ∆x = 1 fm. Dynamical applications, oscillations and collisions, are in general less demanding and can be performed very well with ∆x = 1 fm. This is pleasing as dynamical calculations are usually much more costly than purely static ones.
There are two parameters regulating the static iteration according to Eq. (12), the damping energy e0inv and the step size x0dmp. e0inv should correspond to the depth of the binding potential. The overall step size x0dmp can be of order of one if e0inv is well chosen. Nuclear binding is very similar    all over the chart of nuclei. This allows to develop one safe choice for nearly all cases. We recommend e0inv≈ 100 MeV together with x0dmp≈ 1/2, reducing the latter slightly if convergence problems appear. Of course, a few percent in iteration speed my be gained by fine-tuning these parameters for a given case, but this is not worth the effort unless large scale surveys for a given class of nuclei and forces are planned. The time stepping using the exponential propagator has the two parameters, step size dt and order mxp= m of the Taylor expansion (18) of the exponential. Intuitively, one expects that small dt and large mxp improve the quality of the step. An efficient stepping scheme, however, looks for the largest dt and smallest mxp which still provide acceptable and stable results. It is hard to give general rules as good working values for the parameters depend on all details of the actual calculation: gridding, nuclei involved, excitation energy, and kind of excitation. Figure 4 demonstrates the dependence of a typical dynamical evolution on these time-stepping parameters. We consider a time interval up to 1260 fm/c which is a long time for heavyion collisions and just sufficient for a spectral analysis of oscillations [41]. The excitation is done by a soft sin 2 pulse of finite extension in time. The energy increases during the initial excitation phase, as can be seen from the left panel in the figure. After the external pulse is over, energy conservation holds, which is nicely seen at plotting resolution in the left panel. Normalization should be conserved at all times. Both conservation laws serve as tests for the time step. Norm is conserved up to at least six digits for all cases and times shown in Figure 4. The energy is more critical. The right panels show the energy in a small window around the final energy after the excitation phase is over. The right lower panel shows a variation of the Taylor order m for fixed time step. The most prominent effect is the sudden turn to catastrophic failure for m = 6. In fact, propagation by approximate exponential evolution explodes sooner or later in all cases. The art is to extend the stable interval by a proper choice of the stepping parameters. It is plausible that the cases with m > 6 maintain stability longer because the exponential is better approximated. It is surprising that m = 4 is also stable over the whole time interval. There seem to be subtle cancellations of error going on. Considering the stable signals, we see very little differences between the cases. One may generally be happy with low m. It is mainly stability demands which could call for larger m. Note that this is not a generic result. Stability for a given test case should be checked once in a while and particularly before launching larger surveys.
The right upper panel in Figure 4 shows results for different dt (as we have seen, the m values are not important as long as we achieve stable results). Here we see a clear dependence on the step size. The energies remain constant in the average. But there are energy fluctuations and these depend sensitively on dt. Smaller dt yields smaller fluctuations. As far as one can read off from the figure, the amplitude of the fluctuations shrink ∝dt 2 . It depends on the intended analysis to which level of precision the time evolution should be driven. A value of dt≈ 0.4 fm/c will be acceptable in most cases because the average trend remains far smaller than the fluctuations. Here also it must be emphasized that this is not a generic number. Forces with lower effective mass (SV-bas has m * /m = 0.9) are more demanding and usually require smaller dt. On the other hand, running propagation without the spin-orbit term allows even larger time steps because the spin-orbit potential is the most critical piece in the mean-field Hamiltonian. The mix of p and r imposes high demands on the numerical representations. We again strongly recommend running a few tests when switching forces or excitation schemes.

Code structure
The code is completely modularized to provide as large a degree of encapsulation as possible in order to ease modification. Most modules read their operating parameters from an associated NAMELIST and have their local initialization routines. In addition, a modern style of programming is used that employs a minimum number of local variables and streamlined array calculations that make the code lines very close to the physical equations being solved.
Here we give a brief overview over the modules and their purpose. There is a comprehensive manual supplied with the electronic version that gives a detzailes description of all modules. The source files containing the modules have the same names, but all in lower case, with an extension of .f90. The higher-level modules are: Main program: It calls initialization routines, sets up the initial wave functions using either harmonic oscillator states or reading wave functions of static Hartree-Fock solutions from given input files (module Fragments). It then calls either statichf from module Static or dynamichf from module Dynamic to run the calculation.
Static: This contains the code for the static iterations statichf and the subroutine sinfo to generate output of the results.
Dynamic: runs the dynamic calculation in dynamichf and generates output in tinfo. Also controls the inclusion of an external excitation implemented in module External.
Densities: calculates the densities and current densities by summing over the single-particle states.
Meanfield: contains the central physics calculation: the computation of the components of the mean field (subroutine skyrme) and the application of the single-particle Hamiltonian to a wave function (subroutine hpsi).
Coulomb: calculation of the Coulomb potential.
Energies: calculation of the total energies and its various contributions.
External: calculation of the action of an external potential or initial collective boost of the wave functions.
Pairs: Implementation of the pairing correlations in the BCS approximation.
Moment: calculation of moments and deformation parameters for the bulk density.
Twobody: attempts to divide up the system into two separated nuclei and to calculate their properties and relative motion.
The lower-level supporting modules are: Params: general parameters used throughout the code.
Forces: defines parameters of the Skyrme force and the pairing interaction and constructs them according to input.
Grids: defines everything associated with the numerical grid and sets it up.
Levels: definition of the single-particle wave functions and elementary operations on them such as derivatives.
Fragments: controls the reading of static wave functions from precomputed data and setting them up in the grid.
Inout: contains the subroutines for I/O of wave functions and densities.
Trivial: defines some very basic operations on wave functions and densities.
Fourier: sets up the transform plans for the FFTW3 package to calculate Fourier transforms of wave functions and densities.
Parallel: This comes in two versions. The source file parallel.f90 contains the routines to handle MPI message passing, while sequential.f90 sets up essentially dummy replacements for sequential or OpenMP mode.
User: contains a sample user initialization code which can be used as a template for more complicated setups.

Parallelization
For both OpenMP [46] and MPI [47] the code can be run in parallel mode. Parallelization for the static mode works in OpenMP but not in MPI: the reason is in the orthogonalization step which is not easily amenable for distributed-memory parallel computation. This should be worked on in the future.
MPI and OpenMP can be used jointly if there are computing nodes with multiple processors.
In both cases parallelization is done over the wave functions. The code applies the time-development operator or the gradient iteration, which use the fixed set of mean-field components, to each wave function, and this can naturally be parallelized. Computing the mean fields and densities by summing up over single-particle wave functions is also easily parallelizable.
The library FFTW3 [48] itself can also run on multiple processors in parallel. This can be used in addition to OpenMP or MPI, but was found to be helpful only in the sequential version of the code.

OpenMP
The application of the subroutine tstep for propagating one wave function for one time step, and of add densities for adding one wave function's contribution to the mean fields is done in parallel loops. The only complicating factor is that the densities, being accumulated in several subsets, must be kept separate using the REDUCTION(+) clause of OpenMP. The summation cannot be done internally in add densities, because for the half time step the wave functions are immediately discarded after adding their contribution to the densities to avoid having to store the full set at half time. Thus only the combined tstep-add densities loop should be parallelized.
The OpenMP program version can be compiled using the appropriate compiler option. A separate Makefile.openmp is provided which just contains the -fopenmp option for the GNU compiler. The number of parallel threads is not set by the code: the user should set the environment variable OMP NUM THREADS to the desired number.

MPI
In principle MPI uses the same technique as OpenMP, parallelizing over wave functions. In this case, however, each node contains only a fraction of the wave functions. This has several consequences: 1. The time-stepping of the wave functions can be done independently on each node, but requires that the densities are broadcast to all nodes after each half or full time step by summing up partial densities from the nodes in subroutine collect densities. 2. The other calculation that uses wave functions directly is that of the single-particle properties. These are calculated on each node for the wave functions present on that node and then collected using subroutine collect sp properties. 3. Only one node must be allowed to produce output. This is regulated by choosing node #0 and setting the flag wflag. 4. The saving of the wave functions is done in the following way: Node zero writes a header file containing the job information on file wffile, then each node writes a separate file wwfile.001, wwfile.002, and so on up to the number of nodes. This avoids having to collect the wave functions on one node. 5. Using these parallel output files as fragment initialization or restart files is handled so flexibly that they can be read into a different nodal configuration or even a sequential run.
The MPI version needs the appropriate compiler and linker calls for the system used. The sequential or OpenMP versions are obtained simply by linking with sequential.f90 instead of parallel.f90, which replaces the MPI calls with a set of dummy routines and sets up the descriptor arrays for the wave function allocation in a trivial way.
The Makefile.mpi shows the procedure; in practice systems differ considerably and the user should look up the compilation commands for his particular system.

Input description
All the input is through NAMELIST and many variables have default values. The NAMELISTs should be in this file in the order in which they are described here, any NAMELIST not used for a particular job may be omitted or left in the input file, in which case it is ignored. The input is from standard input, so if the data are prepared in a file inputdata and the large output listing is to go into output, the code should be run using, e. g., spinfile: time-dependent total, orbital, and spin angularmomentum data as three-dimensional vectors.
extfieldfile: time dependence of expectation value of the external field.

Namelist force force
This defines the Skyrme force to be used. In most cases it should just use two input values: name : the name of the force, referring to the predefined forces in forces.data.
pairing : the type of pairing, at present either NONE for no pairing, VDI for the volume-delta pairing, or DDDI for density-dependent delta pairing. The pairing parameters are included in the force definition. Note that the pairing type must be written in upper case.
There is also the possibility for inputting a user-defined force; this is described in detail with module Forces in the online technical documentation.

Namelist main main
This contains general variables applicable to both static and dynamic mode. They are mostly defined in module Params. tcoul: determines whether the Coulomb field should be included. Default is true.
trestart: if true, restarts the calculation from wffile. Default is false.
tfft: if true, the derivatives of the wave functions, but not of the densities, are done directly through FFT. Otherwise matrix multiplication is used, but with the matrix also obtained from FFT. Default is true. mprint: control for printer output. If mprint is greater than zero, more detailed output is produced every mprint iterations or time steps on standard output.
mplot: if mplot is greater than zero, a printer plot is produced and the densities are dumped every mplot time steps or iterations. Default is 0.
mrest: if greater than zero, a wffile is produced every mrest iteration or time step. Default is 0.
writeselect : selects the output of densities by giving a string of characters choosing them (see subroutine write densities for details. Default is 'r', i. e., only the density is written. write isospin : determines whether the densities should be output isospin-summed (false) or separately for neutrons and protons (true). Default is false.
nof : (number of fragments) selects the initialization. nof=0: initialization from harmonic oscillator, only for the static case; nof<0: user-defined initialization by subroutine init user in module User; nof>0: initialization from fragment data as determined in NAMELIST fragments.

Namelist grid grid
This defines the properties of the numerical grid.
nx, ny, nz : number of grid points in the three Cartesian directions. They must be even numbers.
dx, dy, dz : spacing between grid points in fm. If only dx is given in the input, all three grid spacings become equal. The grid positions are then set up to be symmetric with the coordinate zero centrally between point number nx/2 and nx/2+1.

Namelist
radinx, radiny, radinz: the radius parameters of the harmonic oscillator in the three Cartesian directions, in fm.
e0dmp: the damping parameter. For its use see subroutine setdmc.
The default value is 100 MeV.
x0dmp: parameters controlling the relaxation. The default value is 0.2. In special cases it may be desirable to change this to accelerate convergence.
serr: this parameter is used for a convergence check. If the sum of fluctuations in the single-particle energies, sumflu goes below this value, the calculation stops. A typical value is 1.E-5, but for heavier systems and with pairing this may be too demanding. rsep : termination condition. If the final state in a two-body reaction is also of two-body character, the calculation is terminated as soon as the separation distance exceeds rsep. Units: fm. No default. The purpose of this variable is to prevent the calulation of continuing into meaningless configurations, like crossing of the boundary.

Namelist
texternal : indicates that an external perturbing field is used. In this case the namelist extern must be present. Default: false.

Namelist extern extern
The variables read here describe the external field that is applied to get the nucleus into a collective vibration. Details can be found in the description of module External. It is read only if the parameter texternal read in namelist dynamic is true.
ipulse : the type of pulse applied. For ipulse=0 the wave function is multiplied with a phase factor that produces an initial excitation.For ipulse=1 a Gaussian time dependence is used, for ipulse=2 a cos 2 one. Default: 0. Details are given in Eqs. (9d) and (9e).
isoext : isospin character of the excitation. If this is zero, protons and neutrons are exited in the same way. For a value of 1, they behave oppositely but with a coupling that leaves the center-of-mass invariant. Default: 0.
tau0, taut : time at which the excitation field reaches its maximum, and width of the pulse. No defaults.
omega : if this is nonzero, the time-dependence of the external field gets an additional cosine factor with frequency omega.
radext, widext : radius and width of a Woods-Saxon-type cutoff factor in radius for the external field. Defaults: 100 fm and 1 fm, which practically implies no damping. Definition in Eq. (9b).
amplq0 : amplitude for quadrupole excitation of the Q 20 type. Defined as usual with respect to the z-axis.

Namelist fragments fragments
The variables in this namelist control fragment initialization for the case of nof>0. Most quantities are dimensioned for the fragments and we indicate this by index "i" in the following.

Namelist user user
This namelist is read only if needed for user initialization (see module User). Its contents depend on the specific user initialization and the only thing to be said here is that it should appear last in the input file. Since the namelist is defined and used only in module User, its name can also be changed arbitrarily, of course.

Output description
The code produces a number of output files containing various pieces of information. The bulky observables, such as densities or currents, are selectively output at certain time steps into special binary output files nnnnnn.tdd), where nnnnnn indicates the iteration or time step number. These files can then be used for further analysis or converted to be used as input in visualization codes. Examples of this are found among the utility codes provided.
The complete set of wave functions is saved at regular intervals of mrest iterations or time steps. Because this leads to large storage requirements, only the last such file in a run is kept. It can be used for restarting the calculation or for inputting fragment wave functions for initializing another calculation.
In MPI mode the wave functions are distributed over several files, each containing only those present on a specific processor. An additional header file contains the remaining information and can be used to read the wave functions even on a different processor configuration.
Aside from these binary files there are a number of text files. The *.res files contain one line for each time step or iteration where output is triggered according to the value of mprint. There is an explanatory header line in these that has a leading '#', so that is treated as a comment by gnuplotthus gnuplot can be used immediately to plot the behavior of any one column of numbers, i. e., the dependence of a physical quantity on iteration or time. In addition more complicated output is printed on standard output, which can be redirected into a file using shell redirection.
The names of the output files can be adjusted using input variables as listed in Section 5.1, so that the files are here denoted by the default names given there. The *.res files are relatively small, so that no mechanism was implemented to suppress them.

File conver.res
This is produced in sinfo only in the static calculation and its purpose is to give a quick impression of the convergence behavior. The numbers given in each line are the iteration count, the total energy in MeV, the relative change in energy from one iteration to the next, the average uncertainties in the singleparticle energies efluct1 and efluct2, the root-mean-square radius in fm and finally the deformation parameters β and γ (see Section 2.6.1). The latter give an impression as to where the nuclear shape is ending up.
For the judging of convergence, the efluct values are more important than the change in total energy, since the energy can remain constant while the wave functions still change considerably.

File monopoles.res
At present this file is generated in subroutine moment shortprint, but only in the dynamic mode. It contains the time, the neutron, proton, and total root-mean-square radii, and the difference of neutron minus proton root-mean-square radii.

File dipoles.res
This file is produced both in the static and dynamic calculations in subroutines sinfo and tinfo. It contains the iteration or time step number followed by the three components of the center of mass vector R and those of the difference of proton minus neutron center-of-mass vectors R (T =1) , both in fm, for the definition see Eq. (21a). The first of these is useful as a check to see whether the center of mass drifts off during the calculation, while the second vector may be useful to look at proton vs. neutron vibrations.

File quadrupoles.res
This is also generated in moment shortprint and thus only in dynamic mode. It contains the Cartesian quadrupole moments (21b) for neutrons, protons, and the full mass distribution followed by the expectation values of x 2 , y 2 , and z 2 for neutrons and protons, all in fm 2 .

File energies.res
This is the important monitoring file for the dynamic calculation. It is written in subroutine tinfo and each line contains the simulation time, the number of neutrons and protons in the system (these should be constant, so this is a stability check), the total energy (again, this should be conserved), the total kinetic energy, and finally the collective energy ecoll separately for neutrons and protons. Units for the energies are all MeV.

Standard output
This contains all the additional information that in most cases is not needed directly for further processing in, e. g., graphics programs. If it should be found necessary to utilize some data from this file, it is in most cases easy to use grep or a scripting language like PERL or PYTHON to extract the necessary data. Of course the code can also be modified to produce additional output files.
The initial part of the output essentially echoes all the data from the NAMELISTs in tabular form, to enable checking the correctness of input data. In the case of fragment initialization this is more involved and is discussed below for the dynamic case, since it is not so common for static calculations.
In general the layout of the information is compact with sequences of "*" characters to provide separation between input groups as some guidance for the eye.

Static calculation
The code first prints the current iteration number. Iteration "0" refers to the state before iterations are started, for the later iteration numbers, the information refers to the end of the iteration.
The an overview of the various energy contributions is printed: the first part is similar to what in in file conver.res, while a second list shows the energies calculated from the density functional and split up for the various contributions.
Next there is a simple printer plot of the density distribution in the (x − z)-plane. This is often quite helpful, since it shows what is going on in the calculation without the need to start a graphics program, which requires converting the data first.
Next there is a listing of single-particle states. For each state this shows its parity, occupation probability wocc (which is called v 2 , as it is interesting mostly in the pairing case), the energy fluctuations sp efluct1 and sp efluct2, the norm, the kinetic and total energies of the state, and finally the expectation values of the three components of the orbital and spin angular momentum, respectively.
Finally a summary of some integrated quantities is given, separately for neutrons, protons, and all nucleons: the particle number, the root-mean-square radius, quadrupole moment, and the average of the coordinates squared, followed by the centerof-mass components.
Then iterations continue and only one line is printed for each, as this may be quite slow and it is important to be able to check progress while the code is running. After mprint iterations the detailed information is repeated.

Dynamic calculation
After echoing the parameters for the dynamic calculation, the fragment definitions are given and all the resulting information is printed: the computed boost values in case of twobody initialization, the properties of the single-particle states read in, including which index in the fragment file is transferred to which index in the total set of wave functions.
In case of an external field, the input data is also echoed in a detailed form.
The time stepping starts and detailed output is produced every mprint steps at the end of the time step. Much of it is similar to the the static case, so only the differences are pointed out.
Because in the dynamic case the situation can have a general three-dimensional character, the full information on the quadrupole tensor (21b) is printed, separately for the neutron, proton, and total mass distributions. The three eigenvalues and associated normalized eigenvectors are given, followed by Cartesian and polar deformation parameters a 0 , a 2 , β, and γ, as defined in Section 2.6.1.
The separation of the two fragments and its time derivative is printed next and repeated every time step, as examining these quantities is meaningful only with more frequent sampling.
The energy information is shortened by omitting the quantities not of interest in the dynamic case. After the printer plot the results of the two-body analysis are given (for the meaning of the various quantities see module Twobody in the online documentation. The Section on "collision kinematics" shows the mass, charge, position, and kinetic energy of the two fragments. It should be kept in mind that the two-body analysis is only valid if the reaction plane is the (x − z)-plane and the results printed may not be useful if the physical situation is not of two-body nature but the code does not recognize that.
The single-particle property list and the integrated quantities are as in the static case, but the energy fluctuations are omitted.
The next time step is then indicated and the fragment separation data are printed for every time step until after mprint steps the full output recurs.

Utilities
A set of short programs is designed to help with further processing of output from the code. The currently available set is described here.
Most of these routines contain a loop to input a file name from the terminal. If it is desired to do this in a loop over a set of files, a simple trick can be used: generate a list of file names using, e. g., ls -1 *.tdd > list and then execute the program with "list" as input, e. g., ./fileinfo < list For Tdhf2Silo a script convert is provided that handles this (see below). It can easily be adapted to the other utilities and must be stored in the same directory as the executable utility program itself in order for the dirname command to work properly.

Fileinfo
This is a short program to print information about binary files generated by Sky3D. It takes the name of either a *.tdd or a wave function file as input and prints out essentially all the information contained in the header. It can be compiled simply by executing gfortran -o fileinfo fileinfo.f90

Inertia
This is intended as an example of an analysis code reading *.tdd files and doing some computation, which can be used as a model for doing similar things. It illustrates looking for the desired field in the file and taking into account whether it is stored as a total density or isospin-separated. Being given a filename as input, it reads the density and calculates the inertia tensor to print all its 9 components. Compile it using gfortran -o Inertia Inertia.f90

Cuts
This utility reads the density from a file nnnnnn.tdd file and produces output files named nnnnnnrxy.tdd, nnnnnnrxz.tdd, and nnnnnnryz.tdd, which contains two-dimensional cuts through the system in the (x, y), (x, z), and (y, z) plane, respectively. The cuts are evaluated at the origin for the third coordinate by averaging the two neighboring planes.
These data files are written in such a format that they can be read by gnuplot for use in its commands for 2-dimensional plotting.
This program is intended again as a template that can be modified for other applications.

Overlap
This is a code to calculate the overlap of two Slater determinants. Given the names of two wave function files (which must contain compatible data: dimensions, force, etc.) it reads the wave functions, generates the matrices of overlaps between one set and the other separately for neutrons and protons, and then calculates the determinant of each, which is the overlap between the two Slater determinants. It prints some summary information: distance between the centers of mass of each set, minimum and maximum diagonal elements, maximum absolute value of off-diagonal elements, and finally the overlaps for protons and neutrons as well as their product.
This code uses subroutines from LINPACK (stored at NETLIB.ORG), which are included in the file det.f with appropriate copyright. It can be compiled using gfortran -o overlap overlap.f90 det.f.

Tdhf2Silo
This program is quite complicated. It reads a set of *.tdd files and converts them into Silo files. Silo is a library for handling scientific datasets developed at Lawrence Livermore National Laboratory (https://wci.llnl.gov/codes/silo/index.h This is the most appropriate library to use in conjunction with the LLNL graphics visualization tool VisIt (https://wci.llnl.gov/ which was found to be highly suitable for plotting Sky3D results and producing movies. The conversion code is quite flexible in that it decides what to produce for the different field types: isospin-summed or not, vector or scalar. They are given appropriate names for Silo with suffixes p and n for protons and neutrons, and x, y, z for the vector components. In the case of vector fields a variable containing the vector definition is also written so that the field can be plotted immediately as a vector field in VisIt.
If the user wants another dataset handling method, the code should be readily adaptable to other libraries. VisIt itself has many ways of importing data, but of course there also alternative 3D visualization systems.

Compilation and linking
To produce executable files the code comes with several Makefiles. The standard Makefile produces a sequential code sky3d.seq, the file Makefile.openmp a parallel code using OpenMP, while Makefile.mpi should produce an MPI distributed system code.
The Makefiles are written for the gfortran compiler and the commands and options must be adapted if other compilers are used. The user may also have to modify the library names and execution of the code under MPI will require consulting the local documentation or system administrator.
No attempt was made to select the compiler and linker options optimally for speed, since experience has shown that optimization at the cutting edge is highly time-dependent. Thus users should do some speed tests before embarking on major calculations.

External libraries needed
The LAPACK library is used in the code to supply the routines ZHBEVD and DSYEV. LAPACK should be installed in most scientific computing centers; if not, the files can be obtained from www.netlib.org and just be added as additional source files to the code. Note that a complete set of routines called by these two subroutines must be downloaded.
The other external routine library that is used is FFTW3. Again, it will be preinstalled in most systems. If not, there are two possibilities: 1. Download the source code from www.fftw.org and compile the library yourself. In our experience this worked smoothly. The generated library can be installed in a system library directory or kept in a user account. In the latter case the use of -lfftw3 in the makefiles does not work anymore and the full path name of the library file must be given. 2. Replace it by another FFT routine. This requires quite a bit of work: FFTW organizes its calls around "plans", which describe a set of operations to be done on the threedimensional arrays. In init fft quite sophisticated plans are set up to, for example, transform in the y-direction for all x-and z-values.This means that all calls to subroutines beginning with dfftw have to be examined and possibly replaced by loops over one-dimensional FFT transforms. This should be relatively straightforward, but there are two more important points to consider: 1) normalization differs between FFT codes. For FFTW transformation followed by inverse transformation multiplies the original data by nx*ny*nz and this factor is taken into account in several places. 2) For the non-periodic case the Fourier transform in the Coulomb solver uses doubled dimensions in all three directions. Some FFT codes have an initialization that sets up the transformation factors depending on the dimension; in such cases the initialization may have to be repeated.

Running with OPENMP
The OPENMP version can be compiled using the file Makefile.openmp, which produces an executable sky3d.omp. The main difference to the sequential makefile is the addition of an openmp compiler option. Since this depends on the compiler used, it may have to be modified. For gfortran the option is -fopenmp, while for Intel Fortran it is simply -openmp.
For controlling the running of the code the user should set the environment variable OMP NUM THREADS to the number of parallel threads to be used (usually the number of processors). In addition it may be necessary to set OMP STACKSIZE. The two parallel loops in dynamichf need to store all the density fields in parallel, and the second loop adds ps4 to that. Taking into account vector fields and isospin, a total of 24 threedimensional COMPLEX(8) fields need to be stored, amounting to nx*ny*nz*24*16 bytes, which is the stack size needed.

Running under MPI
The situation for MPI is a bit more complex than for OPENMP, so that the file Makefile.mpi will almost certainly have to be modified. One crucial difference to the other makefiles is that the module Parallel is now generated from the source file parallel.f90. In addition the compilation commands have to be adapted; something like mpif90 will be needed but is installation dependent. In addition a command like mpirun will be needed for execution; the user is advised to consult local documentation.

Required input
Here it is just summarized what input is needed for a static or dynamic calculation. A full description can be found with the documentation for the NAMELISTs.

Static calculation
The NAMELISTs needed are, in that order: files, force, main, grid, static. In addition, if initialization is from fragments, fragments, and for user initialization possibly user (only if the user initialization requires input).

Dynamic calculation
The NAMELISTs needed are, in that order: files, force, main, grid, dynamic, For external field excitation extern, and in all cases fragments.

Test cases
To allow checking the proper behavior of the code, we provide three test cases exercising different functions: a static calculation for the ground state of 16 O, a dynamic calculation using an external excitation to stimulate a giant resonance in 16 O, and finally a sample deep-inelastic collision of two 16 O nuclei. The test cases directory contains a more detailed description of these cases and what to look for principally in the results.

Caveats concerning the code
The user should be aware of the limitations of the code in various respects but also note some less straightforward procedures for improving accuracy in certain cases.

Static calculations
Since the main use of the code is expected to be in timedependent calculations, the static part is less highly developed. The omission of all symmetry restrictions, while very useful for innovative applications with time-dependence, can cause some problems in the static case.
1. The spin is not aligned along a fixed direction. Since for even-even nuclei Kramers degeneracy operates, two degenerate levels will mix in an uncontrolled way to produce an arbitrary spin alignment. This can be remedied by diagonalizing the spin operators in such a two-level subspace. 2. The center of mass may move away from the origin during the iterations. This is typically a very small effect and will be corrected when the wave functions are placed at a given position in the dynamic initialization. The calculation of observables, however, should always use coordinates relative to the real center of mass. 3. In very heavy nuclei sometimes even a rotation was seen, as the reoccupation of high-lying levels can change the geometric orientation. If this is a serious problem, constraints should be introduced or another code used for the static calculation.
4. The harmonic oscillator initialization can be quite deficient for heavier nuclei. It is planned to develop a Nilssonmodel alternative; meanwhile in case of problems the use of wave functions from axial or symmetry-restricting codes could be implemented by generating a wave function file from such results. Since this will involve interpolation, a number of static iterations should then be performed using the present code to improve stability.

Dynamic calculations
Here it is important what the required accuracy will be. Most exploratory calculations will not pose high accuracy demands. There are several possible ways in which accuracy can be improved if needed: 1. The initial configuration may be improved. If the fragment nuclei are not situated at grid points, one should run a number of static iterations with each fragment situated alone in the new grid. 2. If the fragments are deformed, the energy estimate from the Rutherford trajectory will not be reliable; in this case it is recommended to run a number of dynamic iterations, observe the change in relative distance, and then correct the boost energies to match the correct velocity of relative motion.

Modifying the code
Since the code was developed with a view for easy modification, in this Section we give some advice on how to add new things to it and how to run simulations.

Modifying the Skyrme force
The code comes with a quite large database of Skyrme force parametrizations together with appropriate pairing parameters built in. This is certainly useful to avoid mistakes in the input by having to indicate only the name of the force. Still there will be a need to add new forces and even forces with a different density functional to it, for which we suggest three different approaches.

Direct parameter input
If a specific parametrization is not expected to be a permanent addition to the code, for example if one or several parameters are varied to study the sensitivity of the results to specific parts of the density functional, the best way is to use the facility for giving the parameters directly in the input. This is triggered by using some force name that is not in the database, in which case the routine expects all parameters to be given in NAMELIST forces.

Expanding the database
At present the database of forces is contained in the file forces.data as a long initialization statement. This has the advantage of readability and avoids having to make a database file accessible in every directory used for code applications.
So a new force which will be used more permanently can be added simply by adding the appropriate lines in forces.data and increasing the number assigned to nforce accordingly. There is a slight danger that the total length of the list will exceed the 255 lines allowed by the Fortran standard; in this case either remove some outdated forces, remove the separation lines of all stars, or if there is still a problem, convert the initialization to DATA statements initializing a smaller number of forces in each case.

Adding new physics to the density functional
This can of course require a lot more modifications to the code. Generally speaking, such a new term will appear not only in the density functional, but also leads to new contributions in the single-particle Hamiltonian and may require new types of densities and currents. In addition, it will involve new force parameters. So in general the following steps will be needed: 1. Parameter definition: The new parameter can be added to the general Force type definition in forces.f90. This is the most logical and systematic way, but we recommend it only if the new physics is to be there permanently, since in this case the whole database has to be updated to give a default value -probably zero -to the new parameter for each existing force. The alternative is to leave this parameter as a separate entity, which can still be a module variable in Forces and be input using the same namelist. Derived parameters should also be calculated here. 2. New densities and currents: everything that can be defined directly as a sum over the occupied single-particle wave functions should be defined and calculated in module Densities. It should be a module variable and allocated during initialization similar to the densities already defined. Then the contributions of the different derivatives of the wave function can be accumulated separately as is done for the existing contributions. 3. Calculation of mean-field components: subroutine skyrme in module Meanfield is the place where the fields appearing in the single-particle Hamiltonian are calculated. The difference to module Densities is that wave functions are not involved directly, but only combinations and derivatives of the densities and currents need to be evaluated. The new fields should be defined and allocated and then calculated in subroutine skyrme. The handling of derivatives again can be imitated based on the existing terms. 4. Single-particle Hamiltonian: The additional contribution to the single-particle Hamiltonian must be calculated in subroutine hpsi of module Meanfield. Code has to added to calculate how the additional terms act on the input wave function pinn and the result has to be added to the output wave function pout. Again, for efficiency the spatial derivatives can be handled in separate loops. 5. Contribution to the energy: subroutine integ energy in module Energies must include an additional contribution of the new term, which in this case means simply computing the expression for the energy functional.
Probably it will be useful to add this up in some new variable (defined as a module variable), so that the contribution of the new term can be printed out together with the other contributions in subroutines sinfo of module Static or subroutine tinfo of module Dynamic, respectively. 6. Output of densities: it may be desirable to write out the new densities, currents, or mean-field contributions into the *.tdd files. To do that, subroutine write densities in module Inout must be modified. A new letter to use for writeselect must be defined and selected in the SELECT CASE statement, and then depending on whether it is a scalar or a vector field, write one density or write vec density is called with a descriptive name given to the field. The complication that occurs here is that these writing routines assume isospin-dependent fields and output either isospin-separate or isospin-summed fields depending on the value of write isospin. If the field has no isospin dependence, the statements used to output wcoul should be imitated.

Analyzing the results in new ways
The code provides quite a large number of physical observables in its output files. For new applications it may be necessary to look at additional ones. There are essentially two ways to implement this: (a) If the new observable depends only on density and mean-field components, the easiest way is to use the *.tdd files, where if necessary more fields can be output by modifying the subroutine write densities.
As an example for reading the *.tdd files, we provide a code calculating the tensor of inertia among the utilities. (b) It becomes more complicated if the wave functions have to be used. Here one possibility is to add a call to a user-written routine at the beginning of the static or dynamic loops, which then can use the array psi in any way desired. This may not be a good option if the calculations are lengthy and the analysis routine may have to be modified several times. In this case it is better to generate a new file name for the wffile each time write wavefunctions is called, similar to the way it is done in write densities, and to do the analysis by reading the wave function files.