Developing New Analysis Programs

Developing a New Program

The ANALYSE package contains a set of library modules and include files that can be linked together to make use of structures in the MCTDH program, e.g. the reading of output files from an MCTDH calculation, and allocation of memory for MCTDH entities. Of course, an ANALYSE program may not need to use these modules.

In order to ease the development of an ANALYSE program there is a template file, template.F . This is a bare-bones program, which should help to set up the automatic memory allocation and use the interface modules. Various choices that need to be made are noted by a comment line:
C!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
In the following, this program is described to highlight what is available in the package, and to illustrate the standard structure. The basic idea is that one (or more) MCTDH output file(s) are read, analysed, and output to a file, or to screen.

Program name

The new program must have a unique name, i.e. program template must be replaced with program newname, and the line commented in. This is necessary as the all_compile script recognises a program by this line.

Include files

There are two include files, global.inc and aglobal.inc that include (through other include files) the variables and logicals required in a program. For reading and writing channels.inc iis required. Other include files may be required for a specific program (e.g. flux).

Set default values

The first stage of a program is to set various default values to the variables and flags. This is done in default for entities related to the MCTDH program, and in adefault for those specific to the ANALYSE programs. If new variables are introduced in a program, and need initialising, this is best done here. This is especially true if the new entity is in a common block in an include file.

Declare standard resources required

The next two sections declare what MCTDH entities are to be used. In the first section the relevant flags are set to .true. if information from a DVR or OPER file are to be read. The second part contains logicals defining what MCTDH entities are to be stored in memory and used in the program, e.g. a wavefunction or set of grid coordinates. How these flags are set will affect the pointer and memory allocation (see below).

The variable ncomp is the number of comparison data sets that are required. If only one set (e.g. a single wavefunction) is used, ncomp=0. If two or more sets are to be compared, this variable must be set to the appropriate value. See Comparison data sets below.

Input information

A subroutine is called to read in and evaluate any input parameters, e.g. file names, options etc. The name of the subroutine should be changed from input to a unique name, or else checking programs such as mutil become confused by the various "input" subroutines in the different programs. Two filenames, filein(1:laein) and fileout(1:laeout) for the input and output files are in the include files.

Open files

The input and output files are opened, and if gnuplot information is requested it is written to the output file.

Read and check system information

The information section of the input file is then read to obtain the system data, e.g. number of degrees of freedom, basis, time intervals etc. A version number is also read which is stored in an array element filever(unit) corresponding to the input file unit. This is to allow the format of the file to be known, and it may be neccesary to use this information in an ANALYSE program if the format is version dependent.

Next the system information from the DVR and OPER files are read, if requested. A check that these files match the system in the input file is then made in subroutine chksyserr. This checking procedure is governed by the flags chkdvr etc.

Memeory Allocation

Before calling the subroutine which will do all the work, the memory required must be allocated. The program has 6 memory blocks available for different data types: mc (complex*16), mr (real*8), mi (integer), ml (logical), ms (complex*8), mf (real*4). The sizes of these blocks must be calculated, and arrays are allocated to these using pointers. A C subroutine is then used to dynamically allocate this memory for the analysis.

Standard arrays requested by setting the appropriate logicals are allocated in zeigausw and memausw. For example, if the wavefunction will be read then lpsi=.true.. Information about the wavefunction (number of single-particle functions, primitive bases, pointers for where the single-particle functions are in the wavefunction array, etc.) has been read from the appropriate psi file. In zeigausw, depending on how the wavefunction has been saved, additional arrays may be required (e.g. a single precision complex array). In memausw, space for the wavefunction is then set up in mc, starting at mcpsi, and space for the complex*8 array spsi is set up im ms, starting at mspsi.

After these standard arrays have been allocated, any program specific arrays can be allocated, e.g. further workspace. If the operator is to be read from the OPER file and used, the subroutine memhpsi must be called after this. NB This routine expects that the work arrays are the last entities in the memory blocks. If this is not the case allocated memory will be overwritten. This might seem unneccesarily complex, but is a result of the MCTDH program structure.

Calling the analysis subroutine

The final step is the memory allocation through the C routine caller, which passes memory to the routine receiver. In this routine, the analysis subroutine can be called, selecting the arrays neccesary from the memory blocks. For example, if a wavefunction is required

call subanalyse(mc(mcpsi))


subroutine subanalyse(psi)

complex*16 psi(dgldim)

Again, the name of the subroutine should be changed from subanalyse to something appropriate.

Reading the DVR or OPER file

Now that the memory has been allocated, it is possible to read the DVR file for particular information by changing the dvrdata flags to .true.. For example, if the grid coordinates are required, and lort = .true. has been set, then the array real*8 ort(ortdim) has been allocated in the mr memory block. Set dvrdata(1)=.true., and then

          call rddvr(ort,rdum,rdum,rdum,zdum,
+                    zdum,idum,idum,idum,rdum,chkdvr)

will return the coordinates in the ort array. The pointer zort(f) points the first element of the grid for the degree of freedom f.

If the operator is required, this can be obtained from the subroutine rdoper, replacing rdum by the real*8 array hops(hopsdim), which has been previously allocated.

Comparison data sets

To enable the MCTDH files to be read with all the various options, and to allow new options to be easily built in, it is neccesary that the various reading routine use the information in the common blocks. Otherwise the argument lists must be changed in all programs. This however makes it difficult to use the same routine for different entities, e.g. if different wavefunctions are to be read. For this reason comparison data sets have been set up.

The information in the common blocks is referred to as the "standard information", e.g. wavefunction array length dgldim, pointers to single-particle functions zetf etc. If more than one wavefunction is to be read the various data sets can be stored in the comparison data sets, and copied to the standard data before use.

To use two data sets, for example, the variable ncomp must be set to 2. After reading the relevant system information this can be stored in a comparison set using the savecX routines. For example

      read(ipsi1) filever(ipsi1)
      chkdvr=1
      chkgrd=1
      chkpsi=1
      chkdat=1
      call rdpsiinfo(ipsi1,chkdvr,chkgrd,chkpsi,chkdat)
      call chksyserr(filename,lerr)
      if (lerr) go to 999

C-----------------------------------------------------------------------
C save comparison set data as set 2
C-----------------------------------------------------------------------
      call savecsys(2)
      call savecpsi(2)

reads the information from the psi file allocated to channel ipsi1, and stroes it as the second data set. A further wavefunction can then be read into set 1. Before reading a wavefunction, this data must be copied to the standard information, e.g.

         call getcpsi(2)
         call rdpsi(ipsi1,psi1,spsi,jindx1)

This comparison data sets have the same names as the standard information, with a c put at the beginning of the array name, and a final index denoting the data set. Thus ctinit(2) is the initial time in the second data set. These comparison sets can thus be directly addressed in a program.

Interface Modules

In this section the modules are described. In order to use them, various flags may have to be set, which are then passed via arguments or common blocks to the routines. There are also dependencies between the routines. ADEFAULT must be called before any modules are used. It may also be neccesary to allocate memory using ZEIGAUSW and MEMAUSW before calling a module. Other useful modules are in the MCTDH program, e.g. keyunits, iofile

ADEFAULT:
Default values are set for the variables used by the modules.

AUSWUTIL:
Contains a number of useful routines. In particular, the routines to copy information between the standard and comparison data sets are here.

MEMAUSW:
Allocates memory for arrays often needed in analyse programs. ZEIGAUSW must be called beforehand. If the logicals listed in the following table are set to true, the memory described is allocated.

logical	Memory allocated
lpsi	wavefunction array
lpsi1	second wavefunction array
lpsigrd	wavefunction arrays in grid representation
ldmat	density matrices
ldmat1	second set of density matrices
lort	array for grid points
lgpop	arrays for grid populations
ltrafo	arrays to transform DVR/FBR.

RDFILES:
This module contains the routines to read the MCTDH output files.
RDPSI returns the wavefunction in MCTDH form.
RDPSIGRID returns the wavefunction in its grid reprentation.
RDGRIDPOP reads the gridpop file.

ZEIGAUSW:
Sets up pointers and array dimensions needed to use standard arrays.