Code Description

Basic Description

The MCTDH program consists of a set of modular routines, connected by include files containing common blocks, and by read-write files. It is important for an understanding of the program structure to realise that a calculation can be divided into four levels, each of which depends on the information from the previous levels:

Input information is read from an input file by the eingabe routines and loaded into common blocks listed in the include files. This information is then passed to the main parts of the program. There are four main parts: rundvr, runoper, runinwf, and runprop. Which parts of the program are to be run depends on the input parameters in the input file, and as to whether the read-write files need to be generated.

Read-Write Files

The structure of the program, and the flow of information between the various parts, and the auxilliary programs POTFIT and ANALYSE, is shown in the figure MCTDH program Structure. The square blocks represent parts of the code, while ellipses represent the read-write files. For each part, referred to as runX, the memory is dynamically allocated. First sizes and pointers for the various arrays are calculated in a routine zeigX.f, or read from the read-write files. The pointers for the arrays are then set up in memX.f. Finally the program part is accessed through a C-routine callX.

The callX routines provide a divide. Information in the include files is above this divide, where the memory for the data arrays has been allocated. This information can be passed through to the runX program parts, but arrays using the dynamically allocated memory cannot pass back through the barrier, and hence the read-write files pass this information below the barrier.

The read-write files all have the following structure:

It is important to note that the first three sections can all be read independently of the data. This allows the information about the data to be retreived before the memory is allocated, as the data can only be read after memory allocation.

The version number is to allow further development. The system definition enables the program to check that the various files match the system of interest, e.g. the primitive basis is correct. The data information is specific to the data stored in the file, e.g. array sizes and pointers.

The system definition is in a set of blocks, as different files need different information. The first block, known as dvrdef, defines the primitive basis DVR, e.g. the number of degrees of freedom, type of basis etc. See wrdvrdef.F for the full definition. The next block, grddef, defines the grid used for the MCTDH calculation, e.g. whether any degrees of freedom are combined together, and how an electronic basis will be treated etc. See wrgrddef.F for the full definition. In the third block, psidef, any further parameters needed to specify the wavefunction are defined, e.g. the number of single-particle functions. See wrpsidef.F for the full definition.

The dvr file contains the results of the DVR basis generated.

The oper file contains the Hamiltonian operator in the primitive basis defined in the dvr file.

The restart file contains the wavefunction.

Include Files

A number of small include files are required to pass the common blocks around the program. These also have a hierachical structure, e.g. the arrays in griddat.inc needed the dimension defined in maxdim.inc. The following list follows this hierachy in the sense that files lower in the list may depend on those higher in the list.

global.inc: Basic information needed for error handeling, lengths of strings etc.

Files containing dimensions needed for arrays in include files:

maxdim.inc: Maximum number of degrees of freedom
maxsta.inc: Maximum number of electronic states
maxkoe.inc: Maximum number of operator terms

Files containing definitions needed by various parts of the program:

timing.inc: Information needed for timing routines
griddat.inc: Information needed to define grid
psidef.inc: Information needed to define the wavefunction
pesrf.inc: Information needed to define potential energy surfaces
operdef.inc: Information needed to define operator
hpsi.inc: Information needed to perform H operation on wavefunction.
daten.inc: Information needed to read the output files.

Files containing information local to a part of the program:

runoper.inc: Information needed by runoper
runinwf.inc: Information needed by runinwf
propwf.inc: Information needed by propwf

In order to use all these files with the minimum of effort, there is an include file for each part of the program which lists all the files required there, e.g. gendvr.inc includes maxdim.inc and griddat.inc, while genoper.inc also includes operdef.inc, pessrf.inc, runoper.inc etc.

Source Code Files

The source code of the MCTDH program is split into several directories, one for each part. This makes the code easier to edit and faster to compile. In comparison to earlier versions some files have been split into smaller files for the abovementioned reasons and some new files have been added. The following list shows all files in alphabetical order (ignoring the library files that are described in the Libraries Description) together with a brief description. To get some detailed information about the subroutines in a specific file click on the file name.

Note that the tables below are currently rather incomplete.

The files in the mctdh directory are:

default.F: Non-local variables and parameters are initialized here.
mdefault.F: MCTDH variables and parameters are initialized here.
eingabe.F: Reads the input data from the input file.
eingutil.F: Contains utility subroutines to open the input and output files and to read the input data.
einrun.F: Reads the Run-Section from the input file.
einsbas.F: Reads the SPF-Basis-Section from the input file.
mctdh.F: MCTDH main program.
rdnatpot.F: Reads natural potential data.
timing.F: Analyses the time needed in the various sections of the program.

The gendvr directory contains files for generating a DVR/FBR representation of the Hamiltonian and single-particle functions:

einpbas.F: Reads the Primitive-Basis-Section from the input file.
init.F: Driver routine for setting up the DVR matrices
init1.F: Initialises the arrays trafo, d1mat and d2mat for a DVR basis set.

The files in the genoper directory serve to build up the Hamiltonian:

einoper.F: Reads the Operator-Section from the input file.
heingabe.F: Reads operator file (cf. Hamiltonian Documentation).
heinutil.F: Auxiliary subroutines to read operator file.
zeigoper.F: Calculates pointer structure and allocates the memory arrays to the standard variable arrays.

In the geninwf directory those files are stored that create an initial wavepacket:

adiab.F: adiabatic correction of the initial state (H+H₂ only)
adiabwkb.F: (diabatic) WKB energy correction
getpsi.F: Contains the subroutines employed to read an old wavefunction.
transmf.F: Computes diabatic translational mean-field.
zeigdvr.F: Calculates pointer structure for rundvr.F (Build of DVR).
zeiginwf.F: Calculates pointer structure for runinwf.F (Build of initial WF).
zeigoper.F: Calculates pointer structure for runoper.F (Build of operator module).

The propwf directory contains files needed for the wavefunction propagation (or the Hamiltonian diagonalisation):

autocorr.F: Contains all subroutines related to the output of the autocorrelation function.
crocofly.F: Contains all subroutines related to the on-the-fly calculation and the output of the cross-correlation function.
daten.F: Contains subroutines needed for data output.
diagham.F: Contains subroutines needed for diagonalising the Hamiltonian.
heval.F: Contains subroutines to calculate the action of the Hamiltonian.
integrat.F: Contains all routines needed for the integration that are MCTDH specific and thus not stored in the various integrator library modules.
matrizen.F: Calculates the time derivative of the MCTDH coefficients. For Hamiltonian terms involving only a single state the mean field matrices are calculated too.
memprop.F: Calculates pointer structure for runprop.
mfields.F: Calculates mean-fields.
output.F: Controls the output of data.
prop.F: Controls the propagation of the wavefunction.
prop1.F: Contains auxiliary subroutines to propagate the wavefunction.
trafos.F: Contains subroutines to transform the wavefunction.
zeigprop.F: Calculates pointer structure and allocates the memory arrays to the standard variable arrays.

The files compile and Makefile are needed to compile the code. For further information see Compiling the Programs.

The file operator.op (where operator stands for the name of the desired operator), specifies the Hamiltonian. For a description see Hamiltonians.

Program Organisation

Variables and Arrays

The program actually starts in the subroutine receiver in mctdh.F. This routine gets six arrays which are dynamically allocated in caller in fort_alloc.c, and whose sizes are determined in zeiger.F. The arrays psi and dtpsi contain wavefunctions, the others - mc, mr, mi, and ml - are general complex, real, integer and logical arrays, and are dimensioned by the parameters maxmemc, maxmemr, maxmemi and maxmeml, respectively. These arrays provide the space for "sub"-arrays used in subroutines, and their structure is given by pointers which are determined in zeiger.F.

The following example shall illustrate the general structure:

The wavefunction is stored in the array psi. As given by the MCTDH-ansatz this includes a (usually) large A-vector and several single-particle functions (SPFs). A common question may be: Where the hell are the single-particle functions for the mode m and the state s? The answer is at psi(zetf(m,s)), i.e. zetf(m,s) points to the beginning of this sub-array. The length of a single SPF is given by subdim(m) and there are dim(m,s) SPFs for one combination (m,s), so the length of our sub-array is dim(m,s)*subdim(m).

We see that the position of a specific array is determined by the corresponding pointer and the corresponding "sub"-array. By this way, the program uses a kind of pseudo-dynamical allocation. The pointers are generally defined and explained in the files global.inc and main.inc, and obtain their values in zeiger.F. Constants are defined in global.inc.

The internal unit system of the program is a.u.

Loop Counters

Standard indices are used for loop counters as follows:

Counter	Description
`m`	counter over modes
`e`	counter over single-particle functions
`f`	counter over degrees of freedom
`k`	counter over Hamiltonian coefficients
`b`	counter over MCTDH coefficients
`s`	counter over states
`g`	counter over grid points
`p`	counter over packets or over parameters
`n`	counter over `ndim`
`v`	counter over `vdim`
`d`	counter over members of psi arrays

This is just to make the code easier to read. For example, psi(b) is then obviously one of the wavefunction coefficients.

Algorithms

In the following some technical details of the various algorithms are described.

CDVR:

The CDVR method serves to employ a potential energy surface without representing it in the product form normally required in the MCTDH scheme. (Description of the CDVR algorithm).

Multi-packet:

In a multi-packet calculation several initial wavepackets are generated and simultaneously propagated. Technically this is done using the multi-set algorithm, i.e. the various packets are considered as nuclear wavefunctions propagated on a set of electronic states, with all these states being described by the same Hamiltonian, and the couplings between them being set to zero. (This mirroring of the Hamiltonian is done automatically by the program.) For example, in the case of the pyrazine model with s = 2 physical electronic states, one then has ps = 2p MCTDH states or sets, where p is the number of packets involved. Note that the current implementation requires equal numbers of single-particle functions for the various packets.

Version number
System definition
Data information
Data