Basic use of the MCTDH package

(Need a general introduction about how the program works including: choice of number of gaussians; what the database is and how to create/maintain it; nos of processors that should be used; definition of final propagation time and propagation time step; use of initial momentum; and 'explain why we would ever need to start from the lower state. Or can we assume that the user should understand this already (and know, for example, how many gaussians is appropriate?))

Creating the Directories

First of all, create a folder in which to save the calculations in your $HOME and $WORK directories. Create a directory in $HOME, and an identical directory in $WORK that has the same path. Keep pathlengths as short as possible (if the total path becomes too long the dynamics code will not run). For simplicity, in this case, we are going to assume that this directory is named buta (for butadiene).

    mkdir buta

Create a subdirectory called dd_data in the new directory in $WORK. Thus you should start with:

   /home/$USER/buta
   /work/$USER/buta

and

   /work/$USER/buta/dd_data

Important Note: it is essential to have the same path in both directories, $HOME and $WORK. Do not change the path-names within the HOME/WORK directories. If you re-name or move any subdirectory in your, for example HOME directory, you have to be consistent and make the same changes into the $WORK directory.

Gaussian Output to be Used as Input for the Dynamics Program

Populate the new directory /work/$USER/buta/dd_data with four files which must be named: coin.log; start.fchk; start.log; and template (no extension). These files are required in order to run any dynamics simulation, they are the basis that the script dd_generator (see the following section) will use to create the necessary files for running MCTDH.

coin.log is the logfile resulting from the optimisation of the conical intersection. NB The MCTDH program will rotate the conical intersection before starting the dynamics. To do this it must be provided with three angles which describe the appropriate rotation. There is an Excel spreadsheet available to calculate these angles (see the section entitled "Rotating the Conical Intersection"). NB the final conical intersection log file must only contain one geometry (otherwise the dynamics program will not know which geometry to use).

Media:coin.log

(Insert a link here to section "Rotating the Conical Intersection")

It is important that the structure of the starting material (in start.fchk and start.log) can be transformed (exactly, including atom labels) into the conical intersection geometry (in coin.log). The results of dynamics calculations will be easier to interpret if the path between the two structures is the shortest possible path and if that path does not involve a high energy barrier. Certain problems can arise if the molecule is highly symmetric. It is easiest to understand the potential difficulties by considering some examples.

Butadiene belongs to the C2h point group. There are four vibrational modes associated with C2h: Ag; Bg; Au; and Bu. Ag is the slope. Travelling from the starting structure to the conical intersection along Ag, the PES for each of the other three modes has a mirror image in the plane of the Ag direction. Thus the wavepacket can go either to the right or to the left (positive or negative) as the symmetry of each of the three other modes is broken. These pathways are equivalent. However, there are consequently (2x2x2 = 2³ =) 8 identical (by symmetry) conical intersections that could be described in coin.log. Obviously coin.log contains only one of these. If, in a dynamics simulation, the wavepacket actually travels to one of the other 7 identical conical intersections, the diabatic surfaces may not be properly described in that region. If an initial momentum is given, the direction (positive or negative) along each coordinate must be specified: this helps to direct the wavepacket towards the desired conical interaction.

Additionally, the starting structure must be able to be transformed exactly into the conical intersection structure via a continuous path, keeping the atomic labels consistent. If the conical intersection geometry has been located through deformation of the starting material and following the PES downhill this is not likely to be a problem (because the optimisation will, in most cases, have followed the shortest path to reach the conical intersection). However, if the conical intersection has been drawn separately, care must be taken to ensure that none of the atom labels have been permuted. For example, for the butadiene molecule shown below, chemical intuition suggests that the first conical intersection should be closer to the starting material on the PES (via a disrotatory motion).

The second conical intersection shown (which is identical by symmetry) is also reached by a disrotatory transformation but it is clear that H10 and H7 have further to travel.

The second conical intersection shown ought to be further away from the starting material on the PES. To model the decay of cis-butadiene through a conical intersection it makes more sense to use the first conical intersection shown above as the reference conical intersection because this conical intersection should be reached more quickly.

start.fchk is the formatted checkpoint file from the high precision vibrational modes of the ground state with state averaged orbitals (Part A, 6)

Media:start.fchk

start.log is the log file from the high precision vibrational modes of the ground state with state averaged orbitals (Part A, 6)

Media:start.log

template is a template gaussian input file to be used by the dynamics program (For example, use the input file from the high precision vibrational modes of the ground state with state averaged orbitals (Part A, 6)). The file template will be the template for the creation of GAUSSIAN input files. It must contain all the mutual keywords of the GAUSSIAN calculations which will be run on the fly by MCTDH.

Media:template.noext

The Generator

The "generator" script uses the Gaussian output files from 2 along with other information supplied by the user to create input files for the dynamics program.

Navigate to /home/$USER/buta and run the generator by typing:

    dd_generator

for the HPC Cluster, or

    dd_generator_MAC

for a MAC OS X system

Typing dd_generator will run the script automatically. After installing the MCTDH package this command has been added to your list of executable commands.

    ********************************************************************************
    ********************    -------  DD GENERATOR  --------    *****************
    ********************************************************************************

The generator will ask a series of questions and will use the answers provided to prepare the input files for the dynamics program. The questions it asks are listed below along with answers for a trial run on butadiene.

List of questions

Question Answer for this tutorial

Which version of GAUSSIAN should we use? (gdv|g03) g03

How many shared-memory processors should we use? (1|2|4|8)
N.B.: (ncpus=1) serial or (ncpus=2|4|8) parallel (omp) version of MCTDH launched as batch job and launching monoprocessor GAUSSIAN interactive jobs (nprocshared=1 or no specification)
1

What is the name of the molecule? butadiene (the first four letters will be used to name the new files that the generator is creating)

What is the number of atoms? 10

What is the number of nuclear Gaussian functions? 1 (1 means 1 gaussian on state 1 and 1 gaussian on state 2)

'What is the final propagation time (in fs)? 100

What is the propagation time step (in fs)? 0.1

From which electronic state should the wavepacket start? (2|1)' 2

Will there be an initial momentum given to the wavepacket? (y|n) n

Do you want to add a reference label to the name of the case? (y|n) n

What will be the status of the database? (rdwr|rd|wr|none) none

>> Manual rotation (using Excel spreadsheet) << Please type the values of the three Euler angles (deg) 180.0000
0.0000
90.0000

Choose rd for read, wr for write etc.
In this case we are not going to use the database, but if it was used (rd or rdwr is specified) the script would ask the user to specify the threshold below which the database values will be used (%). In this case the value none is chosen so this question is not asked.

Please type the value of the maximum-difference criterion (in %): 1 (if 1 is specified, the database value will be used if the geometry is within 1% of the new geometry being calculated)

Files created by the generator.

At the end the script dd_generator tells you that has read the files start.fchk, start.log, and coin.log from the folder /work/dm107/dyn/dd_data/ and written the file start.chk.

Read formatted file /work/dm107/dyn_90dev/dd_data/start.fchk
 Write checkpoint file /work/dm107/dyn_90dev/dd_data/start.chk
   about to read parameters...
   ... parameters just read
   about to read file /work/dm107/dyn_90dev/dd_data/start.log                   
                                                               ...
   ... file /work/dm107/dyn_90dev/dd_data/start.log                             
                                                      read successfully
   about to read file /work/dm107/dyn_90dev/dd_data/coin.log                    
                                                               ...
   ... file /work/dm107/dyn_90dev/dd_data/coin.log                              
                                                      read successfully

and that the extraction of data and the creation of a summary file were done successfully

Data extracted successfully
The summary of this file generation is in:
-rw-r--r--  1 dm107 hpc-users 3070 Dec  2 17:06 /home/dm107/dyn_90dev/but1dd1o.txt
 
Do not forget you still can add a momentum, change the integrator, change convergence criteria for GAUSSIAN and MCTDH...
 
Do not forget to use a re-orientated conical intersection in /work/dm107/dyn_90dev/dd_data/coin.log
 
Do not forget to check the order of the atoms of the conical intersection geometry in /work/dm107/dyn_90dev/dd_data/coin.log
 
Now go to your jobscript file and give values to the memory and walltime.
 
Good luck!

The dd_generator script has generated four files in the same directory (/home/$USER/buta):

but1dd1o_trial.inp
but1dd1o_trial.job
but1dd1o_trial.txt
but1dd_wr.op

(Link to a more detailed explanation of this files)

The script has also created some files and directories inside /work/$USER/buta but where the MCTDH and GAUSSIAN calculations will be stored:

1) Directories:

butadd1o_trial
butadd1o_trial/dd_data

2) Files (Inside butadd1o_trial/dd_data):

refdb.dat
start.chk
template

3) Files inside /work/$USER/buta/dd_data

backward.dat forward.dat

Modifications to be made to the files before running the dynamics.

In the .job file in /home/$USER/buta specify the memory and walltime if you are in the HPC Cluster.

If you wish to add a momentum, specify it in the .inp file in /home/$USER/buta.

In the case of running a development version of GAUSSIAN (gdvg03 or gdvg01), copy the start.chk file into the folder /work/$USER/dyn_90dev/but1dd1o/dd_data. The unformated checkpoint file form start.fchk that dd_generator has created is corrupted.

Running the dynamics.

Navigate to /home/$USER/but. Queue the job file. In this case, type:

    qsub butadd1o.job

for MAC OS X users the shell command is already written in butadd1o.job

Analysis of the results.

()

Detailed information about files and directories

The log file

This file contains information such as the source code version, type of calculation performed, integrator used, numerical parameters, which data files are opened, any error messages, and much more "garbage". The information provided by the log file can be very helpful, in particular when searching for errors. One should always carefully inspect the log file.

The output file

The output contains some standard results:

Time: time in fs. CPU: CPU time in fs. Norm: E-tot: the total energy. This value should be relatively conserved. E-corr: correlated Hamiltonian energy Delta-E: diference between E-tot at time t=0 and the current time t

Note: In a multi-packet run, i.e. when npacket > 1, the total energy and the norm of the wavefunction, as given in the "total" part, are averaged over the packets.

For each state:

pop.: state population E-corr: E-tot

For every mode:

Position expectation value standard deviation Sqrt[<q**2>-**2]

Momentum expectation value

standard deviation Sqrt[<p**2>-

**2] At the end of the propagation, the total CPU-time, host, date and time, path of the name-directory and (if specified in the input) the title of the run are printed.