Resgrp:comp-photo-dyn/mctdh90.56dv/run
Basic use of the MCTDH (DD-vMCG) package
A script called dd_generator uses GAUSSIAN output files and user input to write the input files for mctdh (DD-vMCG).
Notes
The number of gaussian wavepackets (requested by the dd_generator script)
Start with one (to keep calculation times short). When an appropriate momentum, time period, time step and integrator have been established, increase the number of wavepackets until the results are not changed by adding more. Calculation times are proportional to the number of wavepackets. To ensure efficient use of computing resources, always specify an even number of wavepackets if using more than one processor (since the number of processors can only be an even number: currently 2, 4 or 8).
The database
The most time-consuming part of a dynamics run is the GAUSSIAN CASSCF calculations. The results of these calculations can be stored in a database and re-used in future runs. The user must choose whether or not to read and/or write a database. Each calculation will write a separate database so the databases from each run must be combined (and any duplicate entries removed (to save space)). LINK HERE TO SECTION ON MANAGING THE DATABASE
Initial momentum
An initial momentum can be specified to direct the wavepacket to a particular area of the PES (e.g. towards a conical intersection, a transition state or minimum geometry). The momentum is given as a geometry (in frequency mass weighted normal coordinates). The Excel spreadsheet used for rotation of the conical intersection can be used to transform the coordiinates of the desired geometry. If the user requests an initial momentum, dd_generator will search for a file called momentum.dat. If this file is found, dd_generator will ask if the user wishes to use the momentum specified in it. The user can then specify how much energy to give (in eV). Thus, the values in momentum.dat need only be relative values (since they will be multiplied by a scaling factor to provide the total number of eV requested). If there is no file called momentum.dat in the appropriate folder, the user will not be asked to specify a total number of eV to give to the wavepacket(s). A momentum can still be added by manually modifying the .inp file after the dd_generator script has been run (replace the zero values in the fifth column of the section of the .inp file entitled INIT_WF-SECTION). However, it is simpler to create the momentum.dat file mainly because if the user edits the .inp file directly, there is no way of knowing what the total excess energy will be: the values listed in .inp also incorporate magnitude (not just relative values as in momentum.dat).
Memory and number of processors
In template.dat, the memory for each GAUSSIAN CASSCF calculation is specified. In the .job file (created by dd_generator from user input) the total memory is specified. Thus, the memory specified in .job must be at least as great as the memory given in template.dat multiplied by the number of processors.
The GAUSSIAN calculations for each wavepacket will run on a single processor. If the user specifies more than one processor then the calculations for each wavepacket can run simultaneously. To ensure efficient use of computing resources, always specify an even number of wavepackets when using more than one processor (since the number of processors can only be an even number: currently 2, 4 or 8).
Symmetry Considerations
The way the diabatisation is performed in mctdh90.31dv sometimes leads to complications involving symmetry. The diabatic states are defined using the conical intersection (or other geometry) described in coin.log. If the molecule is highly symmetric, there will be more than one conical intersection (all identical by symmetry). If the wavepacket moves close to one of the other conical intersections that is not described in coin.log, the diabatic states may not be properly described in this area of the PES. This of course means that the behaviour of the wavepacket(s) is likely to be described poorly also.
If the conical intersection has lower symmetry than the starting material, the symmetry must necessarily be broken during the dynamics calculation before the conical intersection is reached. The breaking of each mode will occur either due to the accumulation of small numerical inaccuracies or if a momentum is given which forces particular modes to break symmetry. Mctdh90.31dv uses frequency-mass-weighted normal coordinates to describe molecular geometries. Each frequency-mass-weighted coordinate corresponds to a vibration of the molecule (and there are therefore 3N-6 coordinates). Displacements of atoms are described by a (or a combination of) vibration(s). Each vibration can be displaced in one of two directions: positive or negative. The 3N-6 frequency-mass-weighted coordinates can be classified by symmetry. Whether a particular symmetry mode is broken in the positive or negative direction will be equally likely if no initial momentum is given. Thus for each mode that is broken during a simulation the symmetry may break in the positive direction or the negative direction. There will be 2x possible conical intersections corresponding to a particular starting geometry, where x is the number of modes that must be broken going from starting material to conical intersection.
Additionally, care must be taken over atom labeling. It must be possible to transform the starting material into the conical intersection without permuting atom labels. If the conical intersection has been optimised by following the MEP downhill from the starting material this is not likely to be a problem. However, bear this in mind if drawing the conical intersection by hand.
Example: Butadiene
Cis-butadiene belongs to the C2h point group and its 3N-6 frequency-mass-weighted coordinates can be classified as Ag, Bg, Au or Bu (see the C2h character table and GAUSSIAN frequency output). The S1/S0 conical intersection has no symmetry. Thus there are four symmetry modes that must be broken going from starting material to conical intersection. There are sixteen (24) equivalent conical intersections that could be specified in coin.log, corresponding to all possible ways of distorting the starting geometry into the conical intersection geometry. The sixteen conical intersections correspond to: two structures resulting from disrotatory motion in one direction followed by out-of-plane twisting (and their enantiomers); two structures resulting from disrotatory motion in the other direction followed by out-of-plane twisting (and their enantiomers); two structures resulting from conrotatory motion in one direction followed by out-of-plane twisting (and their enantiomers); and two structures resulting from conrotatory motion in the other direction followed by out-of-plane twisting (and their enantiomers). Each of the sixteen conical intersections will have a different arrangement of atomic labels but all can be reached by distorting the starting geometry.
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/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 in 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.dat.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 (DD-vMCG). You may also create a file named "'momentum.dat"' if you wish to give a momentum to the wavepacket(s) (this is optional and may be created later instead).
coin.log is the logfile resulting from the optimisation of the conical intersection. NB The MCTDH (DD-vMCG) 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 structure in coin.log need not necessarily be a true conical intersection: this point is only used for performing the diabatisation. Thus the diabatisation will be exact at this point and will become more approximate the further you are from it (on the PES). If you wish to use a point that is not a conical intersection, simply run a conical intersection optimisation (with only a single step) on the desired geometry. This will create a file with the same layout as a conical intersection optimisation output file (mctdh looks for certain things in this file so if some other type of .log file were used instead, it might not contain the required information). NB the final conical intersection log file must only contain one geometry (otherwise the dynamics program will not know which geometry to use). IS THIS STATEMENT TRUE FOR THIS VERSION??
(Insert a link here to section "Rotating the Conical Intersection")
start.fchk is the formatted checkpoint file from the high precision vibrational modes of the ground state with state averaged orbitals (with atom types rather than numbers) (Part A, 6)
start.log is the log file from the high precision vibrational modes of the ground state with state averaged orbitals (with atom types rather than numbers) (Part A, 6)
template.dat is a template for the creation of GAUSSIAN input files by the MCTDH (DD-vMCG) program. It contains all the mutual keywords of the GAUSSIAN calculations which will be run on-the-fly by MCTDH (DD-vMCG). It also contains other information relating to the GAUSSIAN calculations such as memory requirements.
momentum.dat is a file containing a list of values to be used as vectors to indicate the momentum to be given to a particular (set of) wavepacket(s). The vectors must be listed in the order of the frequency-mass-weighted normal coordinates specified in the Excel spreadsheet. They give the magnitude of the momentum to be given along each of these frequency-mass-weighted normal coordinates. N.B. The signs of the values here are important!
The Generator
- The "generator" script uses the GAUSSIAN output files (above) along with other information supplied by the user to create input files for MCTDH (DD-vMCG).
- 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 (DD-vMCG) 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 asked by dd_generator
Question Answer for this tutorial Explanation 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 The number of nuclear degrees of freedom is 24. Do you want to reduce the dimensionality? (y|n) (N.B.: this can be used also for re-ordering the coordinates)
n What is the number of nuclear Gaussian functions? 1 1 means 1 gaussian on state 1 and 1 gaussian on state 2. This is the "single-set" implementation of MCTDH (DD-vMCG) therefore the gaussian wavepacket on state 1 is the same as the gaussian wavepacket on state 2. 'What is the final propagation time (in fs)? 200 What is the propagation time step (in fs)? 0.1 What is the number of electronic states? (1|2)' 2 What is the highest electronic root?' 2 From which electronic state should the wavepacket start? (2|1)' 2 Will there be an initial momentum given to the wavepacket? (y|n) y Do you want to add a reference label to the name of the case? (y|n) y Please type your text: trial This question is only asked if the answer to the previous question was yes. What will be the status of the database? (rdwr|rd|wr|none) none 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 >> Manual rotation (using Excel spreadsheet) << Please type the values of the three Euler angles (in deg) 176.543006 185.772905
4.538147From the spreadsheet Media:spreadsheet.xls The generator will now tell you that it has read the files start.fchk from the folder /work/dm107/dyn/dd_data/ and written the file start.chk. Read formatted file /work/$USER/buta/dd_data/start.fchk Write checkpoint file /work/$USER/buta/dd_data/start.chk using /apps/gaussian/g03_e01/g03/unfchkDo you want to read the direction of the initial momentum from the existing file in /work/$USER/buta/dd_data? (y|n) (you still have time to change it now, before answering)
y This question is only asked if the answer to the question "Will there be an initial momentum given to the wavepacket?" was yes. Please type the value of the excess kinetic energy (in eV) for calculating the magnitude: 1
Files created by the generator
- The generator will now tell you that it has read the files start.log, coin.log and momentum.dat (if a momentum has been specified) from the folder /work/dm107/dyn/dd_data/:
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
...
*********************************************************
Mass-weighted displacement from start to coin:
old total norm: 22.070387
new total norm: 4.669035
vibrational part: 4.669035
rotational part: 0.000000
*********************************************************
... file /work/dm107/dyn_90dev/dd_data/coin.log
read successfully
about to read file /work/dm107/dyn_90dev/dd_data/momentum.dat
...
*********************************************************
Using the momentum given in momentum.dat would have given
an excess kinetic energy (in eV) of: 15.269
For your information, the momentum has been multiplied
by a factor of: 0.572
*********************************************************
... file /work/dm107/dyn_90dev/dd_data/momentum.dat
read successfully
and that the extraction of data and the creation of a summary file were done successfully. If the Mass-weighted displacement from start to coin appears as NaN it may be because the start.log contains atom numbers rather than atom types. You will need to redo the freq=hpmodes job with atom types.
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 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):
- but1dd1p_trial.inp
- but1dd1p_trial.job
- but1dd1p_trial.txt
- but1dd_none.op
(Link to a more detailed explanation of this files)
- The script has also created some files and directories in /work/$USER/buta where the MCTDH (DD-vMCG) and GAUSSIAN calculations will be stored:
1) Directories:
- butadd1p_trial
- butadd1p_trial/dd_data
2) Files (in butadd1o_trial/dd_data):
- refdb.dat
- start.chk
- template.dat
- forward.dat
- backward.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.
- When running a development version of GAUSSIAN (gdvg03 or gdvg01), copy the start.chk file into the folder /work/$USER/dyn_90dev/but1dd1p_trial/dd_data. The unformatted checkpoint file from start.fchk that dd_generator has created is corrupted (due to a bug in formchk in gdvg01).
Starting from a point other than the Franck-Condon point
When the molecule is highly symmetric and there is more than one equivalent conical intersection (see above) it may be necessary to direct the wavepacket towards the particular conical intersection described in coin.log. The easiest way to do this is to use an initial momentum but it can also be achieved by distorting the geometry of the starting structure slightly towards the conical intersection (i.e. breaking the symmetry in the appropriate direction(s)).
Define the new starting geometry in frequency-mass-weighted normal coordinates (this can be done using the Excel spreadsheet). One option for selecting the new geometry is to use a fraction of the displacement from the FC point to the conical intersection.
After running dd_generator in the normal way, replace the default FC geometry in the .inp file with the new geometry. This geometry (frequency-mass-weighted-coordinates of the centre of the wavepacket) forms the fourth column in the section entitled INIT_WF-SECTION (in .inp).
If the new geometry is very different from the FC geometry, run a single-point state-averaged-CAS calculation (start from RHF, look at the MOs, alter if needed, and check the final active MOs) and call the resulting .chk file start.chk. Replace the start.chk that would have been used ($WORK/buta/but1dd1p_trial/dd_data/start.chk) with this new start.fchk. If the distorted geometry is close enough to the FC geometry it may not be necessary to generate a new start.chk. (When the distortion is large there is a risk of breaking the active space from FC to the new point. Creating a new start.chk is just a way of making sure the starting active space is OK at the new point.) If you are creating a new start.fchk, run the solver on that geometry (relative to the original starting geometry) in the Excel spreadsheet. Create the new start.chk using the geometry in the orange box (so that it is in the same orientation as the molecule in start.log).
N.B. There is no need to replace start.log even when start.fchk is replaced.
Running the dynamics
- Navigate to /home/$USER/buta. Queue the job file. In this case, type:
qsub butadd1p_trial.job
for MAC OS X users the shell command is already written in butadd1p_trial.job
The Database (DB)
The DB is made of four binary files: geo.db pes.db gra.db hes.db They can be used in place of GAUSSIAN calculations to save some computation time when you ask for DBrd or DBrdwr.
Each set of geometry, energy, gradient and Hessian (the latter three quantities having been calculated by GAUSSIAN for the two adiabatic electronic states) will be called a DB record and assigned a record number defining the order of occurrence in the DB. If you run a calculation where the DB is to be read, in the *.res file that you obtain after the calculation is done (file that actually is a redirection of the standard output) you will see the number of records in the DB as it grows. Example (formaldehyde, multi-set 1+1):
53
s= 1 e= 1 t= 0.000 V1= 0.001 V2= 4.533 H11= 0.001 H22=
4.533 H12= 0.000 q= 0.0000 0.0000 0.0000 0.0000
0.0000 0.0000
53
s= 2 e= 1 t= 0.000 V1= 0.001 V2= 4.533 H11= 0.001 H22=
4.533 H12= 0.000 q= 0.0000 0.0000 0.0000 0.0000
0.0000 0.0000
53
s= 1 e= 1 t= 0.500 V1= 0.197 V2= 4.701 H11= 0.211 H22=
4.687 H12= 0.251 q= -0.6258 -0.6146 0.6136 0.6049
-0.5520 -0.0616
53
s= 2 e= 1 t= 0.500 V1= 0.000 V2= 4.533 H11= 0.000 H22=
4.533 H12= -0.001 q= 0.0009 0.0009 0.0093 -0.0245
-0.0066 0.0001
53
ERROR in subroutine Distdb :
Redundant records in Data Base: 18 19
s= 1 e= 1 t= 1.000 V1= 0.469 V2= 4.455 H11= 0.479 H22=
4.445 H12= 0.200 q= -0.4640 -0.5641 1.2594 -0.5871
-1.0864 -0.1120
53
s= 2 e= 1 t= 1.000 V1= 0.001 V2= 4.473 H11= 0.001 H22=
4.473 H12= -0.002 q= 0.0003 0.0004 0.0404 -0.0944
-0.0292 0.0001
53
s= 1 e= 1 t= 1.500 V1= 0.113 V2= 4.543 H11= 0.120 H22=
4.536 H12= 0.170 q= -0.3994 -0.4377 0.5444 0.2057
-0.5149 -0.0638
53
s= 2 e= 1 t= 1.500 V1= 0.013 V2= 4.390 H11= 0.013 H22=
4.390 H12= -0.005 q= 0.0009 0.0010 0.0897 -0.2104
-0.0634 0.0002
53
s= 1 e= 1 t= 2.000 V1= 0.068 V2= 4.431 H11= 0.071 H22=
4.428 H12= 0.110 q= -0.2751 -0.3274 0.4796 -0.0812
-0.3880 -0.0605
54
s= 2 e= 1 t= 2.000 V1= 0.019 V2= 4.367 H11= 0.019 H22=
4.367 H12= -0.009 q= 0.0012 0.0013 0.1563 -0.3658
-0.1058 0.0003
54
s= 1 e= 1 t= 2.500 V1= 0.064 V2= 4.273 H11= 0.065 H22=
4.272 H12= -0.064 q= 0.2559 0.0914 0.1691 -0.7218
-0.1052 -0.0595
55
s= 2 e= 1 t= 2.500 V1= 0.070 V2= 4.258 H11= 0.070 H22=
4.257 H12= -0.013 q= 0.0013 0.0015 0.2384 -0.5572
-0.1528 0.0004
Here, we started with a DB of 53 records. It is large enough until 2 fs, where an additional record is produced, and so on. The error message saying here that records 18 and 19 are identical will not prevent the code from working. It is just here to remind you that the DB is not clean in that it has redundant records. This is to be prevented as the DB files might become large quite rapidly and all records will be read at each update time (dd_update = tout in the integrator section, tout = xxx fs in the run section of the *.inp file) for each Gaussian function.
Incidentally, when you want to analyse results written to a *.res file, you need to look for lines starting with "s=" (other lines are various messages and numbers of DB records). The criterion for using the DB is to be seen as a maximum-allowed relative difference
So, it is up to you to clean the DB when you see this, and to keep this in mind when you concatenate DBs that have branched in their past (see more details below).
Also, a given DB depends on the choice of a reference geometry in Cartesian coordinates because geometries written to the DB files will always be re-orientated with respect to this particular geometry. The back and forth transformations are performed by MCTDH (DD-vMCG) so you don't need not worry about this. This geometry is written to refdb.dat. It is the starting geometry that comes from start.log (the only difference being that it is mass-centred). As long as you do not change the starting geometry or the nature of the electronic states, you can thus keep growing the DB of your system.
MANAGING THE DATA BASE
A set of DB programs can be found in XXXXX. Their roles are explained below. For now, the sources are in dd_utildev and dd_util.31dv and the executables are in mctdh90dev/bin and mctdh90.31dv/bin, so you can call them from anywhere, as long as DB files are there.
1) countdb.out: returns the number of records of the current DB 2) cleandb.out: "brute-force" cleaning of the current DB. This might be too long to be run interactively when the number of records is too large, because each record needs to be compared to each other record (e.g. a DB with ca. 3800 records will take a couple of minutes to clean interactively). If so, run in batch mode (ask me or Marta if you have a doubt). The cleaned files are: geoclean.db pesclean.db graclean.db hesclean.db You need to rename them without 'clean' for them to be used as a DB. 3) cutdb.out: allows you to cut out the beginning of a DB. This is to be used in the following situation: - You start from a mother DB, DB0, assumed to be clean, with N0 records. - You want to use DB0 in two independent calculations. You replicate DB0 to daughter DBs, DB1 and DB2. - When these calculations are done, DB1 and DB2 have grown to N0+N1 and N0+N2 records, respectively. - You want to concatenate them to create a new mother DB, DB3. - DB1 and DB2 have their N0 first records the same. To avoid redundant records, you need to cut N0 records from DB2 for example. --> go to where DB2 is and use cleandb.out The resulting files are: geocut.db pescut.db gracut.db hescut.db these have N1 records. --> rename DB1, for instance geo.db1 etc., and concatenate the files: cat geo.db1 geoclean.db > geo.db etc. - Your new mother DB, DB3, should now have N0+N1+N2 distinct records. 4) deldb.out: allows you to remove one selected record (for instance if corrupted) from the DB if you know its number. The resulting files are: geodel.db pesdel.db gradel.db hesdel.db You need to rename them without 'del' for them to be used as a DB. 5) seedb.out: allows you extract one selected record from the binary files for reading it (the record is not removed from the DB). The resulting file can be read as a text file (ascii) and is: see.db
Restarting dynamics jobs
The examples below make use of the directory and filenames in the examples above.
A dynamics run can be restarted if the job was terminated before completion (e.g if it ran out of time or failed because of convergence problems) or if the user wishes to continue the job for a longer time period than had originally been specified. The program will continue writing the same output file but will overwrite the .res file. Thus the .res file should be saved under a different name to prevent loss of the data calculated in the original run. The .res is not written until the job stops so this can be done after the restart.
- To increase the final propagation time
Retrieve the .inp file (in $HOME/buta) and change propagation to continuation in the section entitled RUN-SECTION. Then change the final time tfinal to a larger value. (NB this new value should include the time already done).
Re-queue the job.
- To increase the number of scf cycles (the most common reason for convergence failure)
Edit the template.dat file in $WORK/buta/but1dd1p_trial/dd_data. Insert scf=(maxcycle=128) on the keyword line (for example). It is important that you change the template.dat file in the but1dd1p_trial/dd_data folder that was created by dd_generator (rather than the one created by the user in $WORK/buta/dd_data) because the original $WORK/buta/dd_data/template.dat file is only read by dd_generator and used to create the job-specific file in the job folder. (Of course if you wish all future jobs run in $WORK/buta to have increased scf cycles, change the other template.dat file in $WORK/buta/dd_data too.)
Re-queue the job.
- To increase the run-time so that a job will finish
Edit the walltime in the .job file in the $HOME directory and re-queue the job.
Analysis of the results
The most important files are the output and .res files.
Detailed information about files and directories
- output
The output contains some standard results:
Time: time in fs.
CPU: CPU time in s.
Norm:
E-tot: the total energy in eV. This value should be conserved (i.e by the end of the run the E-tot value should not be too different from what it was a time=0).
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.
population : diabatic populations of state 1 and state 2 (in that order) at the current time t
For every mode:
< q >: position expectation value
<dq>: standard deviation Sqrt[< q2 > - < q >2]
< p > : momentum expectation value
<dq>: standard deviation Sqrt[< p2 > - < p >2]
- .res
This file contains the energy of each state (diabatic and adiabatic) at each timestep for each gaussian function:
s is a redundant value in mctdh90.31dv. It is only meaningful in multi-set implementations of mctdh (where there are double the number of gaussian functions than specified: one for each state). (It indicates whether the gaussian function in question is on the lower (1) or upper (2) state. (In the single set implementation, the gaussian functions on each state are identical)).
e is a label given to each gaussian function
t is the current time in fs
V1 is the energy of the lower adiabatic state at time t
V2 is the energy of the upper adiabatic state at time t
H11 is the energy of the lower diabatic state at time t
H22 is the energy of the upper diabatic state at time t
H12
q is the list of frequency mass-weighted normal coordinates that describe the geometry at time t
The end of the file will read something like:
Propagation was successful.
Total time [h:m:s] : 0 : 1 : 21.60
---------------------------------------------------------------------------
------ Host: "cx1-3-2-4" ----------Sat Jan 17 06:50:20 2009
The other most common result is that the calculation aborted because of scf convergence problems in the CASSCF GAUSSIAN calculations. In this case the ouput file will break off after the last complete step and the end of the .res file will read something like:
ERROR in subroutine Getddpes :
CAS convergence failed for gwp on state: 1 1
If the job terminates due to scf convergence problems, it may be possible to restart it after modifying the .inp and/or .job and/or template.dat files.
- log
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. 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.
- dvr
- oper
- ddpeserr
- input
- op.log
- psi
- check
- auto
- ddpes
- restart
- stop
- speed
- timing
- update
Example: butadiene
- Examine the end of the output file to check that the calculation ran to completion.
- Plot the energy of the diabatic and adiabatic states over time (.res file: V1, V2, H11 and H22). Calculate V2-V1 and H22-H11 and plot these also. H12 may also be plotted on the same graph.
The diabatic states cross at 176.2 fs and the adiabatic states cross at 174.5 fs. At a conical intersection, the value V2-V1 will become negative and the H22-H11 will be a minimum (<1eV). Ideally these crossings will occur at the same time. If the diabatisation has worked well, after V2-V1 becomes negative, V2-V1 and H22-H11 will become mirror images (i.e. the diabatic and adiabatic states follow each other).
- Plot the value of the total energy over time (.res file: E-tot). Verify that the Delta-E value is not too large (compare with the total energy (E-tot) at time 0).
- Plot the population of one of the (diabatic) states over time (.res file: population).
At the conical intersection the population transfer is large. In this case the job has been terminated after 200 fs. Usually the calculation should be continued until the populations stabilise (which obviously has not happened in this case).
Making a Movie
Navigate to /work/$USER/buta/but1dd1p_trial (the place where the .res file was written). Type:
dd_movie.out
- The script will ask a series of questions and will use the answers provided to create a file that can be read into Molden or iMol.
List of questions asked by dd_movie.out
Question Answer for this tutorial Explanation What is the number of atoms? 10 What is the actual number of active coordinates? 24 In this case all coordinates were active (so there are 3N-6 where N is the number of atoms) Name of result file to read? but1dd1p_trial.res Ensure that the forward.dat and backward.dat files are located in the current directory in the folder dd_data (they will have been written there automatically) Name of movie file to write? but1dd1p_trial.xyz Choose any name .xyz Which electronic state (s=...)? 1 Choose 1 (this question dates from the multi-set implementation of mctdh where the wavepackets on each state were different) Which Gaussian function (e=...)? 1 In this case we have only one wavepacket The script will now return the number of lines and the initial number of points number of lines: 4684 initial number of points: 4674Maximum number of points you want? 4674 If you wish to view the movie in Molden you must reduce the number of points to less than 100 (so choose 99 for maximum number of points). iMol can cope with any number of frames so choose the number printed by the script: this number will be reduced and the script will tell you how many points it has written.
The .xyz file is ready to be opened in Molden or iMol.
Current Known Bugs
Parallelisation:
It is not possible to run calculations in parallel with version 90.31dv (reading the same template.dat file simultaneously by each processor is causing problems). In 90.35tmp this is ok. In 90.4X it will also be fixed.
Total energy:
1) The total energy is sometimes "artificially" increased by a lot. If so, try a different integrator or play around with the timestep (may need to get larger or smaller). Or try increasing scf convergence to 10-12.
2) If Delta-E suddenly becomes too large, you can get the following message: "ERROR in subroutine WRGAUSSIAN :
Do not use $Swap:...$ metastring in file /.../template.dat"
Often, when getting this message, no more steps are added to the simulation, although the calculation seems to keep running (the integrator keeps trying and failing).
Using the database:
1) Every record generated with mctdh90.31dev and several GWPs is corrupted, so, for now, when launching calculations with more that one wavepacket, use "formdd_none.op" (i.e. do not read nor write the database). [For instance, when using dbrdwr and starting from a clean db, all newly generated records are added to the existing database, even if they are identical to records that were already written in the db.] In 90.35tmp this is ok. In 90.4X it will also be fixed.
2) Records obtained in calculations with one GWP should be OK, but this is still to be confirmed.
"Impossible case":
Error message that forces the calculation to stop. It seems to be associated with moving towards an "unexpected" geometry (e.g. a radical).
State Averaging:
To increase the threshold above which stateaveraging is turned off, increase the default value of 10 (to e.g. 50) in the line ddregions in the OPERATOR-SECTION of the .inp file.
Questions and Modifications to be made to this wiki
- Why we would ever need to start the gwps from the lower state? To study GS rxns; chemiluminesence (eg glowworms).
- In terms of Delta-E, how much is too much? ca1eV but shape of E-tot vs time is important: discontiunities or v flat then steep increase are less reliable than smooth gradual increases.
- Benjamin: Note that Marta has noticed that for very strange geometries, you may receive an error message saying "impossible case", which stops MCTDH. It is too complicated to explain here. I know where this happens in the code, and I should prevent this error to stop the code but I didn't have time to fix this. I'll do it soon.
- Need to edit/correct section on the DB.
- Is does unfchk in gdvg03 create corrupted fchk files? G03 is fine. we don't know about gdvg03 but we don't care because we probablly won't use it.