You can download the manual as pdf or read the html version below.
SHARC: Surface Hopping in the Adiabatic Representation Including Arbitrary Couplings – Manual
Contents
1 Introduction
1.1 Capabilities
1.2 References
1.3 Authors
1.4 Suggestions and Bug Reports
1.5 Notation in this Manual
2 Installation
2.1 How To Obtain
2.2 Terms of Use
2.3 Installation
2.3.1 Libraries
2.3.2 Test Suite
2.3.3 Additional Programs
2.3.4 Quantum Chemistry Programs
3 Execution
3.1 Running a single trajectory
3.1.1 Input files
3.1.2 Running the dynamics code
3.1.3 Output files
3.2 Typical workflow for an ensemble of trajectories
3.2.1 Initial condition generation
3.2.2 Running the dynamics simulations
3.2.3 Analysis of the dynamics results
3.3 Auxilliary Programs and Scripts
3.3.1 Setup
3.3.2 Analysis
3.3.3 Interfaces
4 Input files
4.1 Main input file
4.1.1 General remarks
4.1.2 Input keywords
4.1.3 Detailed Description of the Keywords
4.1.4 Example
4.2 Geometry file
4.3 Velocity file
4.4 Coefficient file
4.5 Laser file
5 Output files
5.1 Log file
5.2 Listing file
5.3 Data file
5.3.1 Specification of the data file
5.4 XYZ file
6 Interfaces
6.1 Interface Specifications
6.1.1 QM.in Specification
6.1.2 QM.out Specification
6.1.3 Further Specifications
6.1.4 Save Directory Specification
6.2 MOLPRO Interface
6.2.1 Interface specific input file: SH2PRO.inp
6.2.2 Template file: MOLPRO.template
6.2.3 Error checking
6.2.4 Things to keep in mind
6.2.5 Molpro input generator: molpro_input.py
6.3 MOLCAS Interface
6.3.1 Interface specific input file: SH2CAS.inp
6.3.2 Template file: MOLCAS.template
6.3.3 Template file generator: molcas_input.py
6.4 COLUMBUS Interface
6.4.1 Template input
6.4.2 Interface specific input file: SH2COL.inp
6.4.3 Template setup
6.5 Analytical PESs Interface
6.5.1 Parametrization
6.5.2 Interfacespecific input file
7 Auxilliary Scripts
7.1 Wigner Distribution Sampling: wigner.py
7.1.1 Usage
7.1.2 Output
7.1.3 Nondefault Masses
7.2 Setup of Initial Calculations: setup_init.py
7.2.1 Usage
7.2.2 Input
7.2.3 Input for MOLPRO
7.2.4 Input for COLUMBUS
7.2.5 Input for MOLCAS
7.2.6 Input for analytical potentials
7.2.7 Input for Run Scripts
7.2.8 Output
7.3 Excitation Selection: excite.py
7.3.1 Usage
7.3.2 Input
7.3.3 Matrix diagonalization
7.3.4 Output
7.3.5 Specification of the initconds.excited file format
7.4 Setup of Trajectories: setup_traj.py
7.4.1 Input
7.4.2 Interfacespecific input
7.4.3 Run script setup
7.4.4 Output
7.5 Laser field generation: laser.x
7.5.1 Usage
7.5.2 Input
7.6 Calculation of Absorption Spectra: spectrum.py
7.6.1 Input
7.6.2 Output
7.7 File transfer: retrieve.sh
7.8 Data Extractor: data_extractor.x
7.8.1 Usage
7.8.2 Output
7.9 Plotting the Extracted Data: make_gnuscript.py
7.10 Calculation of Ensemble Populations: populations.py
7.10.1 Usage
7.10.2 Output
7.11 Obtaining special geometries: crossing.py
7.11.1 Usage
7.11.2 Output
7.12 Internal Coordinates Analysis: geo.py
7.12.1 Input
7.12.2 Options
7.13 Diagonalization Helper: diagonalizer.x
8 Methodology
8.1 Absorption Spectrum
8.2 Active and inactive states
8.3 Damping
8.4 Decoherence
8.5 Excitation Selection
8.6 Internal coordinates definitions
8.7 Kinetic energy adjustments
8.8 Laser fields
8.8.1 Form of the laser field
8.8.2 Envelope functions
8.8.3 Field functions
8.8.4 Chirped pulses
8.8.5 Quadratic chirp without Fourier transform
8.9 Laser interactions
8.9.1 Surface Hopping with laser fields
8.10 Random initial velocities
8.11 Representations
8.11.1 Current state in MCH representation
8.12 Sampling from Wigner Distribution
8.13 Scaling
8.14 Seeding of the RNG
8.15 Selection of gradients and nonadiabatic couplings
8.16 State ordering
8.17 Surface Hopping
8.18 Velocity Verlet
8.19 Wavefunction propagation
8.19.1 Propagation using nonadiabatic couplings
8.19.2 Propagation using overlap matrices
Chapter 1 Introduction
When a molecule is irradiated by light, a number of dynamical processes can take place, in which the molecule redistributes the energy among different electronic and vibrational degrees of freedom. Kasha’s rule [1] states that radiationless transfer from higher excited singlet states to the lowestlying S_{1} excited state is faster than fluorescence (F). This radiationless transfer is called internal conversion (IC) and involves a changes between electronic states of the same multiplicity. If a transition occurs between electronic states of different spin, the process is called intersystem crossing (ISC). A typical ISC process is from a singlet to a triplet state, and once the lowest triplet is populated, phosphorescence (P) can take place. In figure 1.1, radiative (F and P) and radiationless (IC and ISC) processes are summarized in a socalled Jablonski diagram.
The nonradiative IC and ISC process are fundamental concepts which play a decisive role in photochemistry and photobiology. IC processes are present in the excitedstate dynamics of many organic and inorganic molecules, whose applications range from solar energy conversion to drug therapy. Even many, very small molecules, for example O2 and O3, SO2, NO2 and other nitrous oxides, show efficient IC, which has important consequences in atmospheric chemistry and the study of the environment and pollution. IC is also the first step of the biological process of visual perception, where the retinal moiety of rhodopsin absorbs a photon and nonradiatively performs a torsion around one of the double bonds, changing the conformation of the protein and inducing a neural signal. Similarly, protection of the human body from the influence of UV light is achieved through very efficient IC in DNA, proteins and melanins. Ultrafast IC to the electronic ground state allows quickly converting the excitation energy of the UV photons into nuclear kinetic energy, which is spread harmlessly as heat to the environment.
ISC processes are completely forbidden in the frame of the nonrelativistic Schrödinger equation, but they become allowed when including spinorbit couplings, a relativistic effect [2]. Spinorbit coupling depends on the nuclear charge and becomes stronger for heavy atoms, therefore it is typically known as a “heavy atom” effect. However, it has been recently recognized that even for molecules with only first and secondrow atoms, ISC might be relevant and can be competitive in time scales with IC. A small selection of the growing number of molecules where efficient ISC in a subps time scale has been predicted are SO2 [3,[4,[5], benzene [6], aromatic nitrocompounds [7] or DNA nucleobases and derivatives [8,[9,[10,[11,[12].
Theoretical simulations can greatly contribute to understand nonradiative processes by following the nuclear motion on the excitedstate potential energy surfaces (PES) in real time. These simulations are called excitedstate dynamics simulations.
Since the BornOppenheimer approximation is not applicable for this kind of dynamics, nonadiabatic effects need to be incorporated into the simulations.
The principal methodology to tackle excitedstate dynamics simulations is to numerically integrate the timedependent Schrödinger equation, which is usually called full quantum dynamics simulations (QD). Given accurate PESs, QD is able to reach or surpass experimental accuracy. However, the need of an “a priori” knowledge of the full multidimensional PES renders this type of simulations quickly unfeasible for more than few degrees of freedom. Several alternative methodologies are possible to alleviate this problem. One of the most popular ones is to use surface hopping nonadiabatic dynamics.
Surface hopping was originally devised by Tully [13] and greatly improved later by the “fewestswitches criterion”[14] and it has been reviewed extensively since then, see e.g. [15,[16,[17].
In surface hopping, the motion of the excitedstate wavepacket is approximated by the motion of an ensemble of many independent, classical trajectories. Each trajectory is at every instant of time tied to one particular PES, and the nuclear motion is integrated using the gradient of this PES. However, nonadiabatic population transfer can lead to the switching of a trajectory from one PES to another PES. This switching (also called “hopping”, which is the origin of the name “surface hopping”) is based on a stochastic algorithm, taking into account the change of the electronic population from one time step to the next one.
The advantages of the surface hopping methodology and thus its popularity are well summarized in Ref. [15]:
 The method is conceptually simple, since it is based on classical mechanics. The nuclear propagation is based on Newton’s equations and can be performed in cartesian coordinates, avoiding any problems with curved coordinate systems as in QD.
 For the propagation of the trajectories only local information of the PESs is needed. This avoids the calculation of the full, multidimensional PES in advance, which is the main bottleneck of QD methods. In surface hopping dynamics, all degrees of freedom can be included in the simulation. Additionally, all necessary quantities can be calculated ondemand, usually called “onthefly” in this context.
 The independent trajectories can be trivially parallelized.
The strongest of these points of course is the fact that all degrees of freedom can be included easily in the calculations, allowing to describe large systems.
One should note, however, that surface hopping methods in the standard formulation [13,[14] – due to the classical nature of the trajectories – do not allow to treat some purely quantummechanical effects like tunneling, (tunneling for selected degrees of freedom is possible [18]). Additionally, quantum coherence between the electronic states is usually described poorly, because of the independenttrajectory ansatz. This can be treated with some adhoc corrections, e.g., in [19].
In the original surface hopping method, only nonadiabatic couplings are considered, only allowing for population transfer between electronic states of the same multiplicity.
The SHARC methodology is a generalization of standard surface hopping since it allows to include any type of coupling. Beyond nonadiabatic couplings (for IC), spinorbit couplings (for ISC) or interactions of dipole moments with electric fields (to explicitly describe laserinduced processes) can be included.
A number of methodologies for surface hopping including one or the other type of potential couplings have been proposed in references [20,[21,[22,[23,[24,[25,[26], but SHARC can include all types of potential couplings on the same footing.
The SHARC methodology is an extension to standard surface hopping which allows to include these kinds of couplings. The central idea of SHARC is to obtain a fully diagonal Hamiltonian, which is adiabatic with respect to all couplings. The diagonal Hamiltonian is obtained by unitary transformation of the Hamiltonian including all couplings. Surface hopping is conducted on the transformed electronic states.
This has a number of advantages over the standard surface hopping methodology, where no diagonalization is performed:
 Potential couplings (like spinorbit couplings and laserdipole couplings) are usually delocalized. Surface hopping, however, rests on the assumption that the couplings are localized and hence surface hops only occur in the small region where the couplings are large. Within SHARC, by transforming away the potential couplings, additional terms of nonadiabatic (kinetic) couplings arise, which are localized.
 The potential couplings have an influence on the gradients acting on the nuclei. To a good approximation, within SHARC it is possible to include this influence in the dynamics.
 When including spinorbit couplings for states of higher multiplicity, diagonalization solves the problem of rotational invariance of the multiplet components (see [24]).
The SHARC suite of programs is an implementation of the SHARC method. Besides the core dynamics code, it comes with a number of tools aiding in the setup, maintenance and analysis of the trajectories.
1.1 Capabilities
The main features of the SHARC suite are:
 Nonadiabatic dynamics based on the surface hopping methodology able to describe internal conversion and intersystem crossing with any number of states (singlets, doublets, triplets, or higher multiplicities).
 Algorithms for stable wavefunction propagation in the presence of very small or very large couplings.
 Inclusion of interactions with laser fields in the longwavelength limit. The derivatives of the dipole moments can be included in strongfield applications.
 Propagation using either nonadiabatic couplings vectors 〈α[(∂)/(∂R)]β〉, timederivative couplings 〈α[(∂)/(∂t)]β〉 or wavefunction overlaps 〈α(t_{0})β(t)〉 (via the local diabatization procedure [19]).
 Gradients including the effects of spinorbit couplings (with the approximation that the diabatic spinorbit couplings are slowly varying).
 Energydifferencebased partial coupling approximation to speed up calculations [27].
 Energybased decoherence correction [19].
 Flexible interface to quantum chemistry programs. Existing interfaces to MOLPRO 2010 and 2012, MOLCAS 7.8 and 8.0, and to COLUMBUS 7.0 (only if interfaced to MOLCAS). An interface to implement dynamics based on analytical expressions is also available.
 Calculation of Dyson norms for singlephoton ionization spectra through the COLUMBUS interface.
 Suite of auxiliary Python scripts for all steps of the setup procedure and for various analysis tasks.
 Comprehensive tutorial.
1.2 References
The following references should be cited when using the SHARC suite:
 [28] M. Richter, P. Marquetand, J. GonzálezVázquez, I. Sola, L. González:
“SHARC: ab Initio Molecular Dynamics with Surface Hopping in the Adiabatic Representation Including Arbitrary Couplings”,
J. Chem. Theory Comput., 7, 12531258 (2011).  [29] S. Mai, M. Richter, M. Ruckenbauer, M. Oppel, P. Marquetand, L. González:
“SHARC: Surface Hopping Including Arbitrary Couplings – Program Package for NonAdiabatic Dynamics”.
sharcmd.org (2014).  S. Mai, P. Marquetand, L. González:
“A general method to describe intersystem crossing dynamics in trajectory surface hopping” ,
Int. J. Quant. Chem., DOI: 10.1002/qua.24891 (2015).
Further details can be found in the following references:
The theoretical background of SHARC is described in Refs. [28,[30,[31,[32].
Applications of the SHARC code can be found in Refs. [9,[11,[33,[4,[34].
Other features implemented in the SHARC suite are described in the following references:
 Energybased decoherence correction: [19].
 Local diabatization and wavefunction overlap calculation: [35,[27,[36].
 Sampling of initial conditions from a quantummechanical harmonic Wigner distribution: [37,[38].
 Excited state selection for initial condition generation: [39].
 Calculation of ring puckering parameters and their classification: [40,[41].
The quantum chemistry programs to which interfaces with SHARC exist are described in the following sources:
1.3 Authors
The current version of the SHARC suite has been programmed by Sebastian Mai and Martin Richter of the AG González of the Institute of Theoretical Chemistry of the University of Vienna with contributions from Jesús GonzálezVázquez, Philipp Marquetand, Matthias Ruckenbauer, Markus Oppel and Leticia González.
1.4 Suggestions and Bug Reports
Bug reports and suggestions for possible features can be submitted to sharc@univie.ac.at.
1.5 Notation in this Manual
Names of programs
The SHARC suite consists of Fortran90 programs as well as Python and Shell scripts. The executable Fortran90 programs are denoted by the extension .x, the Python scripts have the extension .py and the Shell scripts .sh. Within this manual, all program names are given in bold monospaced font.
Chapter 2 Installation
2.1 How To Obtain
SHARC can be obtained from the SHARC homepage www.sharcmd.org. In the Download section, register with your email adress and affiliation. You will receive a download link to the stated email adress. Clicking on the link in the email will download the archive file containing the SHARC package. Note that the link is active only for 24 h and the number of downloads is limited.
Note that you must accept the Terms of Use given in the following section in order to download SHARC.
2.2 Terms of Use
1. Recitals
The University of Vienna (in following the “University”), whose principal address is at Universitätsring 1, 1010 Vienna, Austria, developed the Software “SHARC: Surface Hopping with Arbitrary Couplings. A program suite to study the excitedstate dynamics of molecules” (in the following the “Software”).
The Software is not to be used for commercial purposes. Our Software provided on the relevant website is solely to be used for research, teaching and demonstration purposes and is free of any charge. The user of the Software is required to have the appropriate knowledge and knowhow to use the Software (in the following the “User”). Any User should seek further guidance and make independent enquiries before relying on the Software.
2. Research, Teaching and Demonstration Purposes, Content and Liability
The Software on this website is based on scientific work and is solely to be used for research, teaching and demonstration purposes. The University, its organs, associated companies, employees and representatives do not warrant, represent or guarantee the functionality of the Software whatsoever.
The use of the Software is at your own risk and your own responsibility. The University, its organs, associated companies, employees and representatives do not assume any liability for the use the Software including its results. The University does not warrant, represent or guarantee any quality or performance of the Software and shall not be responsible for including but not limited to
 any contents of the Software in any kind;
 any intermediate or permanent amendments of the Software of any kind;
 any damages of any kind arising out of or in connection with the use of the Software (including without limitation damages for any consequential loss or loss of business opportunities or projects, or loss of profits) including, but not limited to, data theft or data loss;
 any errors, faults, incompleteness, inaccuracy, inadequacy or unsuitability or incorrectness of the provided Software or its results;
 any personal damages or damages to property due to the use of the Software;
 any damages of computers, any viruses or Trojans or similar which were transferred or conveyed by the Software;
 any defaults or incompleteness of the Software including any loss or damages of any kind due to or arising out of the use of the Software provided on the website of the University;
 any damages caused by third persons or subcontractors.
Hence, the University explicitly reserves the right to amend, supplement, delete or remove the Software, partially or as a whole, at any time and without prior notice.
The University explicitly refrains from and does not assume any liability whatsoever for any products or services offered or advertised by third persons via the provided Software.
3. Acknowledgements
The User shall acknowledge the authors or originators and use of the Software in the publication of any results achieved through use of the Software in the form of including the following references:
 S. Mai, M. Richter, M. Ruckenbauer, M. Oppel, P. Marquetand, L. González: “SHARC: Surface Hopping Including Arbitrary Couplings – Program Package for NonAdiabatic Dynamics”. www.sharcmd.org (2014).
 M. Richter, P. Marquetand, J. GonzálezVázquez, I. Sola, L. González: “SHARC: ab Initio Molecular Dynamics with Surface Hopping in the Adiabatic Representation Including Arbitrary Couplings”. J. Chem. Theory Comput., 7, 12531258 (2011).
The User shall reproduce a copyright notice on every copy of the Software including partial copies and on any accompanying manuals and documentation in the form “Copyright © 2014 University of Vienna. All rights reserved.” Trademark and other proprietary notices must also be reproduced but the User has no longer the right to use the name, arms, trademark, logo or other designation of the University of Vienna.
4. Intellectual Property
The contents of this Software are exclusively governed and construed by Austrian intellectual property laws. Any copying, editing, distribution or use of the contents require the prior consent of the respective authors or originators which is not granted or deemed to be granted without the prior written consent as described above.
If and to the extent the User amends, edits or otherwise supplements the Software with own information or other information of third persons, the User shall be deemed as the new operator of the Software. Consequently, the User expressly assumes any and all liability in connection with the Software.
5. Indemnity
The User is irrevocably obliged to indemnify the University, its employees, its organs, its associated companies and representatives for any damages in particular claims and damage claims including legal cost for
 any use of the Software by the User;
 any violation or infringement of this terms of use partially or in whole;
 any violation of rights of third persons in particular of intellectual property rights or personality rights of publicity or
 any claims of third persons due to the use of the Software.
This obligation of the User to indemnify shall not impede any other provision of this disclaimer described herein and shall remain in force and effect after the User has ceased to use the Software.
If and to the extent third persons assert any claims against the University due to or arising out of the use by any User, such User shall be obliged to indemnify and hold harmless the University for any damages asserted against the University.
6. Data Protection
The User explicitly acknowledges that the University uses services of third party suppliers (e.g. Google Analytics etc.) and consents to the use of its data by the third party supplier. Further, the User consents to the use of the respective IP address of the User or similar by third persons (e.g. for statistical use). The University shall be entitled to provide third persons with such data for statistical purposes. If personal data of Users is provided to or used by third persons, the University shall not be deemed as a sponsor pursuant to the Austrian Data Protection Act 2000.
7. Final Provisions
The Software was developed in Austria and is applicable through a server located in Austria. Any and all interpretation of the content of the Software, any claims with regards to the Software and any disputes in connection with the Software shall be governed by Austrian law excluding any applicable collision provisions. Exclusive place of jurisdiction shall be Vienna, Austria.
2.3 Installation
In order to install and run SHARC under Linux (Windows and OS X are currently not supported), you need the following:
 A Fortran90 compiler (This release is tested against GNU Fortran 4.4.7 and Intel Fortran 15.0).
 The BLAS, LAPACK and FFTW3 libraries.
 Python 2 (This release is tested against Python 2.6.6 and 2.7.3).
The source code of the SHARC suite is distributed as a tar archive file. In order to install it, first extract the content of the archive to a suitable directory:
tar xzvf sharc.tgz
This should create a new directory called sharc/ which contains all the necessary subdirectories and files.
To compile the Fortran90 programs of the SHARC suite, go to the source/ directory.
cd source/
and edit the Makefile by adjusting the F90 variable to point to the Fortran compiler of your choice.
Issuing the command:
make
will compile the source and create all the binaries.
make install
will copy the binary files into the sharc/bin/ directory of the
SHARC distribution, which already contains all the python scripts which
come with SHARC.
In figure 2.1 the directory structure of the complete SHARC directory is shown.
In order to use the SHARC suite, set the environment variable $SHARC to the bin/ directory of the SHARC installation. This ensures that all programs of the SHARC suite find the other executables and all calls are successful. For example, if you have unpacked SHARCinto your home directory,
just set:
export SHARC=~/sharc/bin (for bourne shell users)
or
setenv SHARC $HOME/sharc/bin (for cshell type users)
Note that it is advisable to put this line into your shell’s login
scripts.
2.3.1 Libraries
SHARC requires the BLAS, LAPACK and FFTW3 libraries. During the installation, it might be necessary to alter the LDFLAGS string in the Makefile, depending on where the relevant libraries are located on your system. In this way, it is for example possible to use vendorprovided libraries like the Intel MKL. For more details see the INSTALL file which is included in the SHARC distribution.
2.3.2 Test Suite
After the installation, it is advisable to first execute the test suite of SHARC, which will test the fundamental functionality of SHARC.
Change to an empty directory and execute
$SHARC/tests.py
The interactive script will first verify the Python installation (no message will appear if the Python installation is fine). Subsequently, the script prompts the user to enter which tests should be executed. Their is at least one test for each of the auxiliary scripts and interfaces. Tests whose names start with scripts_ test the functionality of the auxiliary programs in the SHARC suite. Tests whose names start with Analytical_, COLUMBUS_, MOLCAS_ or MOLPRO_ run short trajectories, testing whether the main dynamics code, the interfaces and the quantum chemistry programs work together correctly.
If the installation was successful and Python is installed correctly, the tests named scripts_<NAME> and Analytical_overlap should execute without error.
The test calculations involving the quantum chemistry programs ( COLUMBUS, MOLCAS, MOLPRO) can be used to check that SHARC can correctly call these programs and that they are installed correctly.
If any of the tests show differences between output and reference output, it is advisable to check the respective files. Note that small differences in the output can already occur when using a different version of the quantum chemistry programs, and that these small differences can add up during the simulations.
2.3.3 Additional Programs
For full functionality of the SHARC suite, several additional programs are recommended:
 The Python package NUMPY.
 The GNUPLOT plotting software.
 A program for molecular visualization, able to read files in the xyz file format (e.g. MOLDEN, GABEDIT, MOLEKEL or VMD)
Optimally, within your Python installation the NUMPY package (which provides many numerical methods, e.g., matrix diagonalization) should be available. If NUMPY is not available, the SHARC suite is still functional, and the affected scripts will fall back to use a small Fortran code (frontend for LAPACK) within the SHARC package. Since in the Python scripts no largescale matrix calculations are carried out, there should be no significant performance loss if NUMPY is not available.
GNUPLOT is not strictly necessary, since all output files could be plotted using other plotting programs. However, a number of scripts from the SHARC suite automatically generate GNUPLOT scripts after data processing, allowing to quickly plot the output files.
2.3.4 Quantum Chemistry Programs
Even though SHARC comes with an interface for analytical potentials (and hence can be used without any quantum chemistry program), the main application of SHARC is certainly onthefly ab initio dynamics. Hence, one of the following interfaced quantum chemistry programs is necessary:
 MOLPRO (this release was checked against MOLPRO 2010 and 2012).
 MOLCAS (this release was checked against MOLCAS 7.8 and 8.0).
 COLUMBUS 7, interfaced to MOLCAS, for correlated multireference wavefunctions.
See the relevant sections in chapter 6 for a description of the quantum chemical methods available with each of these programs.
Chapter 3 Execution
The SHARC suite consists of the main dynamics code sharc.x and a number of auxiliary programs, like setup scripts and analyse tools. Additionally, the suite comes with interfaces to suitable quantum chemistry software, e.g. MOLPRO, MOLCAS or COLUMBUS.
In the following, first it is explained how to run a single trajectory by setting up all necessary input for the dynamics code sharc.x manually. Afterwards, the usage of the auxiliary scripts is explained. Detailed infos on the SHARC input files is given in chapter 4 and on the auxiliary scripts in chapter 7. The interfaces are described in chapter 6.
3.1 Running a single trajectory
3.1.1 Input files
SHARC requires a number of input files, which contain the settings for the dynamics simulation (input), the initial geometry (geom), the initial velocity (veloc), the initial coefficients (coeff) and the laser field (laser). Only the first two (input, geom) are mandatory, the others are optional. The necessary files are shown in figure 3.1.
The content of the main input file is explained in detail in section 4.1, the geometry file is specified in section 4.2. The specifications of the velocity, coefficient and laser files are given in sections 4.3, 4.4 and 4.5, respectively.
Additionally, the directory QM/ and the script QM/runQM.sh need to be present, since the onthefly ab initio calculations are implemented through these files. The script QM/runQM.sh is called each time SHARC performs an onthefly calculation of electronic properties (usually by a quantum chemistry program). In order to do so, SHARC first writes the request for the calculation to QM/QM.in, then calls QM/runQM.sh, waits for the script to finish and then reads the requested quantities from QM/QM.out. The script QM/runQM.sh is fully responsible to generate the requested results from the provided input.
In virtually all cases, this task is handled by the SHARCquantum chemistry interfaces (see chapter 6), so that the script QM/runQM.sh has a particularly simple form:
cd QM/
$SHARC/<INTERFACE> QM.in
with the corresponding interface name given. Note that the interfaces in all cases need additional input files, which must be present in QM/. Those input files contain the specifications for the quantum chemistry information, e.g., basis set, active and reference space, memory settings, path to the quantum chemistry program, path to scratch directories; or for SHARC_Analytical.py, the expressions for the analytical potentials. For each interface, the input files are slightly different. See sections 6.2, 6.3, 6.4 or 6.5 for the necessary information.
3.1.2 Running the dynamics code
Given the necessary input files, SHARC can be started by executing
user@host> $SHARC/sharc.x input
Note that besides the input file, at least the geometry file needs to be present (see in the input section for details).
3.1.3 Output files
Figure 3.2 shows the content of a trajectory directory after the execution of SHARC. There will be six new files. These files are output.log, output.lis, output.dat and output.xyz, as well as restart.ctrl and restart.traj.
The file output.log contains mainly a listing of the chosen options and the resulting dynamics settings. At higher print levels, the log file contains also information per timestep (useful for debugging). output.lis contains a table with one line per timestep, giving active states, energies and expectation values. output.dat contains a list of all important matrices and vectors at each timestep. This information can be extracted with data_extractor.x to yield plottable table files. output.xyz contains the geometries of all timesteps (the comments to each geometry give the active state).
For details about the content of the output files, see chapter 5.
The restart files contain the full state of a trajectory and its control variables from the last successful timestep. These files are needed in order to restart a trajectory at this timestep (either because the calculation failed, or in order to extend the simulation time beyond the original maximum simulation time).
3.2 Typical workflow for an ensemble of trajectories
Usually, one is not interested in running only a single trajectory, since a single trajectory cannot reproduce the branching of a wavepacket into different reaction channels. In order to do so, within surface hopping an ensemble of independent trajectories is employed.
When dealing with a (possibly large) ensemble of trajectories, setup and analysis need to be automatized. Hence, the SHARC suite contains a number of scripts fulfilling different tasks in the usual workflow of setting up ensembles of trajectories.
The typical workflow is given schematically in figure 3.3.
3.2.1 Initial condition generation
In the typical workflow, the user will first create a set of suitable initial conditions. In the context of the SHARC package, an initial condition is a set of an initial geometry, initial velocity, initial occupied state and initial wavefunction coefficients.
Many such sets are needed in order to setup physically sound dynamics simulations.
Generation of initial geometries and velocities
Currently, within the SHARC suite, initial geometries and velocities can be generated based on a quantum harmonic oscillator Wigner distribution. The theoretical background is given in section 8.12. The calculation is performed by wigner.py, which is explained in section 7.1.
As given in figure 3.3, wigner.py needs as input the result of a frequency calculation in MOLDEN format. The calculation can be performed by any quantum chemistry program and any method the user sees fit.
wigner.py produces the file initconds, which contains a list of initial conditions ready for further processing.
Generation of initial coefficients and states
In the second preparation step, for each of the sampled initial geometries it has to be decided which excited state should be the initial one. In simple cases, the user may manually choose the initial excited state using excite.py (see 7.3). Alternatively, the selection of initial states can be performed based on the excitation energies and oscillator strengths of the excited states at each initial geometry (this approximately simulates a deltapulse excitation).
The latter option makes it necessary to carry out vertical excitation calculation before the selection of the initial states.
The calculations can be set up with setup_init.py (see section 7.2). This script prepares for each initial condition in the initconds file a directory with the necessary input to perform the calculation. The user should then execute the run script (run.sh) in each of the directories (either manually or through a batch queueing system).
After the vertical excitation calculations are completed, the vertical excitation energies and oscillator strengths of each calculation are collected by excite.py (see 7.3). The same script then performs the selection of the initial electronic state for each initial geometry. The results are written to a new file, initconds.excited. This file contains all information needed to setup the ensemble.
Additionally, spectrum.py (7.6) can calculate absorption spectra based on the initconds.excited file. This may be useful to verify that the level of theory chosen is appropriate (e.g., by comparing to an experimental spectrum), or to choose a suitable excitation window for the determination of the initial state.
3.2.2 Running the dynamics simulations
Based on the initial conditions given in initconds.excited, the input for all trajectories in the ensemble can be setup by setup_traj.py (see section 7.4). The script produces one directory for each trajectory, containing the input files for SHARC and the selected interface.
In order to run a particular trajectory, the user should execute the run script (run.sh) in the directory of the trajectory. Since those calculations can run between minutes and several weeks (depending on the level of theory used and the number of timesteps), it is advisable to submit the run scripts to a batch queueing system.
The progress of the simulations can be monitored most conveniently in the output.lis files. If the calculations are running in some temporary directory, the output files can be copied to the local directory (where they were setup) with the scp wrapper retrieve.sh (see 7.7). This allows to perform ensemble analysis while the trajectories are still running.
If a trajectory fails, the temporary directory where the calculation is running is not deleted. The file README will be created in the trajectory’s directory, giving the time of the failure and the location of the temporary data, so that the case can be investigated.
3.2.3 Analysis of the dynamics results
Each trajectory can be analyzed independently by inspecting the output files (see chapter 5). Most importantly, calling data_extractor.x (7.8) on the output.dat file of a trajectory creates a number of formatted files which can be plotted with the help of make_gnuscript.py (7.9) and GNUPLOT.
The nuclear geometries in output.xyz file can be analyzed in terms of internal coordinates (bond lengths, angles, ring conformations, etc.) using geo.py (7.12).
For the complete ensemble, the script populations.py (section 7.10 can calculate average excitedstate populations. The script crossing.py (7.11) can find and extract notable geometries, e.g., those geometries where a surface hop between two particular states occured.
3.3 Auxilliary Programs and Scripts
The following tables list the auxiliary programs in the SHARC suite. The rightmost column gives the section where the program is documented.
3.3.1 Setup
wigner.py  Creates initial conditions from Wigner distribution.  7.1 
setup_init.py  Sets up initial vertical excitation calculations.  7.2 
excite.py  Generates excited state lists for initial conditions and selects initial states.  7.3 
setup_traj.py  Sets up the dynamics simulations based on the initial conditions.  7.4 
laser.x  Prepares files containing laser fields.  7.5 
molpro_input.py  Prepares MOLPRO input files for common jobs.  6.2.5 
molcas_input.py  Prepares template files for the MOLCAS interface.  6.3.3 
diagonalizer.x  Helper program for excite.py. Only required if NUMPY is not available.  7.13 
3.3.2 Analysis
spectrum.py  Generates absorption spectra from initial conditions files.  7.6 
retrieve.sh  scp wrapper to retrieve dynamics output during the simulation.  7.7 
data_extractor.x  Extracts plottable results from the SHARC output data file.  7.8 
make_gnuscript.py  Creates gnuplot scripts for a given number of states.  7.9 
populations.py  Calculates ensemble populations.  7.10 
crossing.py  Extracts specific geometries from ensembles.  7.11 
geo.py  Calculates internal coordinates from xyz files.  7.12 
3.3.3 Interfaces
SHARC_MOLPRO.py  Allows for the calculation of SOC, gradients, nonadiabatic couplings, time derivatives, overlaps, dipole moments and angular momentum expectation values at the CASSCF level of theory. Symmetry or RASSCF are not supported. Only segmented basis sets are possible.  6.2 
SHARC_MOLCAS.py  Allows for the calculation of SOC, gradients, overlaps and dipole moments at the CASSCF and RASSCF level of theory. Symmetry is not supported.  6.3 
SHARC_COLUMBUS.py  Allows for the calculation of SOC, gradients, overlaps, dipole moments and Dyson norms at the CASSCF, RASSCF, MRCISD and LRTMRAQCC levels of theory. Symmetry is not supported. Only works with the COLUMBUS– MOLCAS interface.  6.4 
SHARC_Analytical.py  Allows to calculate SOC, gradients, overlaps, dipole moments and dipole moment gradients based on analytical expressions of diabatic matrix elements. Works only with cartesian coordinates.  6.5 
Chapter 4 Input files
In this chapter, the format of all SHARC input files are presented. Those are the main input file (here called input), the geometry file, the velocity file, the coefficients file and the laser file. Only the first two are mandatory, the others are optional input files. All input files are ASCII text files.
4.1 Main input file
This section presents the format and all input keywords for the main SHARC input. Note that when using setup_traj.py, full knowledge of the SHARC input keywords is not required.
4.1.1 General remarks
The input file has a relatively flexible structure. With very few exceptions, each single line is independent. An input line starts with a keyword, followed optionally by a number of arguments to this keyword. Example:
stepsize 0.5
Here, stepsize is the keyword, referring to the size of the timesteps for the nuclear motion in the dynamics. 0.5 gives the size of this timestep, in this example 0.5 fs.
A number of keywords have no arguments and act as simple switches (e.g., restart, gradcorrect, grad_select, nac_select, ionization, track_phase, dipole_gradient). Those keywords can be prefixed with no to explicitly deactivate the option (e.g., norestart deactivates restarts).
In each line a trailing comment can be added in the input file, by using the special character #. Everything after # is ignored by the input parser of SHARC. The input file also can contain arbitrary blank lines and lines containing only comments. All input is caseinsensitive.
The input file is read by SHARC by subsequently searching the file for all known keywords. Hence, unknown or misspelled keywords are ignored. Also, the order of the keywords is completely arbitray. Note however, that if a keyword is repeated in the input only the first instance is used by the program.
4.1.2 Input keywords
In Table 4.1, all input keywords for the SHARC input file are listed.
Table 4.1: Input keywords. The first column gives the name of the keyword, the second lists possible arguments and the third line provides an explanation. Defaults are marked like this. $n denotes the nth argument to the keyword.


Keyword  Arguments  Explanation 
printlevel 
integer  Controls the verbosity of the log file. 
$1=0  Log file is empty  
$1=1  + List of internal steps  
$1=2  + Input parsing information  
$1=3  + Some information per timestep  
$1=4  + More information per timestep  
$1=5  + Much information per timestep  
restart  Dynamics is resumed from restart files.  
norestart  Dynamics is initialized from input files.  
norestart takes precedence.  
nstates  list of integers  Number of states per multiplicity. 
$1 (1)  Number of singlet states  
$2 (0)  Number of doublet states  
$3 (0)  Number of triplet states  
$… (0)  Number of states of higher multiplicities  
actstates  list of integers  Number of active states per multiplicity. 
same as nstates  By default, all states are active.  
state  integer, string  Specifies the initial state (no default; SHARC exits if state is missing). 
$1  Initial state.  
$2=MCH  Initial state and coefficients are given in MCH representation.  
$2=diag  Initial state and coefficients are given in diagonal representation.  
coeff  string  Sets the wavefunction coefficients. 
$1=auto  Initial coefficient are determined automatically from initial state.  
$1=external  Initial coefficients are read from file.  
coefffile  quoted string  File containing the initial wavefunction coefficients. 
“coeff”  
geomfile  quoted string  File name containing the initial geometry. 
“geom”  
veloc  string  Sets the initial velocities. 
$1=zero  Initial velocities are zero.  
$1=random $2 float  Initial velocities are determined randomly with $2 eV kinetic energy per atom.  
$1=external  Initial velocities are read from file.  
velocfile  quoted string  File containing the initial velocities. 
“veloc”  
laser  string  Sets the laser field. 
$1=none  No laser field is applied.  
$1=internal  Laser field is calculated at each timestep from internal function.  
$1=external  Laser field for each timestep is read during initialization.  
laserfile  quoted string  File containing the laser field. 
“laser”  
laserwidth  float  Laser bandwidth used to detect induced hops. 
1.0 eV  
stepsize  float  Length of the nuclear dynamics timesteps in fs. 
0.5 fs  
nsubsteps  integer  Number of substeps for the integration of the electronic equation of motion. 
25  
nsteps  integer  Number of simulation steps. 
3  
tmax  float  Total length of the simulation in fs. 
No effect if nsteps is present.  
killafter  float  Terminates the trajectory after $1 fs in the lowest state. If $1 < 0, trajectories are never killed. 
1  
surf  string  Chooses the potential energy surfaces used in surface hopping. 
$1=sharc  Uses the diagonal potentials.  
$1=fish  Uses the MCH potentials.  
coupling  string  Chooses the quantities describing the nonadiabatic couplings. 
$1=ddr  Uses vectorial nonadiabatic couplings 〈ψ_{α}∂/∂Rψ_{β}〉.  
$1=ddt  Uses temporal nonadiabatic couplings 〈ψ_{α}∂/∂tψ_{β}〉.  
$1=overlap  Uses the overlaps 〈ψ_{α}(t_{0})ψ_{β}(t)〉 (Local Diabatization).  
gradcorrect  Includes (E_{α}−E_{β})〈ψ_{α}∂/∂Rψ_{β}〉 in the gradient transformation.  
nogradcorrect  Transforms only the gradients itself.  
ekincorrect  string  Adjustment of the kinetic energy after a surface hop. 
$1=none  Kinetic energy is not adjusted. Jumps are never frustrated.  
$1=parallel_vel  Velocity is rescaled to adjust kinetic energy.  
$1=parallel_nac  Only the velocity component in the direction of 〈ψ_{α}∂/∂Rψ_{β}〉 is rescaled.  
grad_select  Only some gradients are calculated at every timestep.  
grad_all  All gradients are calculated at every timestep (Alias: nograd_select).  
grad_all takes precedence.  
nac_select  Only some 〈ψ_{α}∂/∂Rψ_{β}〉 are calculated at every timestep.  
nac_all  All 〈ψ_{α}∂/∂Rψ_{β}〉 are calculated at every timestep (Alias: nonac_select).  
nac_all takes precedence.  
eselect  float  Parameter for selection of gradients and nonadiabatic couplings (in eV). 
0.5 eV  
select_directly  Selection of gradients and NACs in one call.  
ezero  float  Energy shift for the diagonal elements of the Hamiltonian (in hartree). 
0.0  Is not determined automatically.  
scaling  float  Scaling factor for the Hamiltonian matrix and the gradients. 
1.0  0. < $1  
dampeddyn  float  Scaling factor for the kinetic energy at each timestep. 
1.0  0. ≤ $1 ≤ 1.  
rngseed  integer  Seed for the random number generator. 
10997279  
decoherence  Applies decoherence correction.  
nodecoherence  No decoherence correction.  
nodecoherence takes precedence.  
decoherence_param  float  The value α in the decoherence correction (in Hartree). 
0.1  0. < $1  
ionization  Calculates ionization probabilities onthefly.  
noionization  No ionization probabilities.  
ionization_step  integer  Calculates ionization probabilities every $1 timestep. 
1  By default calculated every timestep (if at all).  
track_phase  Follows the phase of the transformation matrix U.  
notrack_phase  No phase following of U (not recommended, only for debugging).  
dipole_gradient  Include the derivatives of the dipole moments in the gradients.  
nodipole_gradient  Neglect the derivatives of the dipole moments. 
4.1.3 Detailed Description of the Keywords
Printlevel
The printlevel keyword controls the verbosity of the log file. The data output file (output.dat) and the listing file (output.lis) are not affected by the printlevel. The printlevels are described in section 5.1.
Restart
There are two keywords controlling trajectory restarting. The keyword restart enables restarting, while norestart disables restart. If both keywords are present, norestart takes precedence. The default is no restart.
When restarting, all control variables are read from the restart file instead of the input file. The only exceptions are nsteps and tmax. In this way, a trajectory which ran for the full simulation time can easily be restarted to extend the simulation time.
Number of States and Active States
The keyword nstates controls how many states are taken into account in the dynamics. The keyword arguments specify the number of singlet, doublet, triplet, etc. states. There is no hardcoded maximum multiplicity in the SHARC code, however, some interfaces may restrict the maximum multiplicity.
Using the actstates keyword, the dynamics can be restricted to some lowest states in each multiplicity. For each multiplicity, the number of active states must not be larger than the number of states. All couplings between the active states and the frozen states are deleted. These couplings include offdiagonal elements in the H^{MCH} matrix, in the overlap matrix and in the matrix containing the nonadiabatic couplings. Freezing states can be useful if transient absorption spectra are to be calculated without increasing computational cost due to the large number of states.
Note that the initial state must not be frozen.
Initial State
The initial state can be given either in MCH or diagonal representation. The keyword state is followed by an integer specifying the initial state and either the string mch or diag. For the MCH representations, states are enumerated according to the canonical state ordering, see 8.16. The diagonal states are ordered according to energy. Note that the initial state must be active.
If the initial state is given in the MCH basis, determination of the initial diagonal state is carried out after the initial QM calculation.
Initial Coefficients
The initial coefficients can be read from a file. The default filename is coeff, but the filename can be given with the keyword coefffile. Note that the filename has to be enclosed in single or double quotes. The file must contain the real and imaginary part of the initial coefficients, one line per state with no blank lines inbetween. These coefficients are interpreted to be in the same representation as the initial state, i.e. the state keyword influences the initial coefficients. For details on the file format, see section 4.4.
Alternatively, the initial coefficients can be determined automatically from the initial state, using coeff auto in the input file. If the initial state is given in the diagonal representation as i, the initial coefficients are c^{diag}_{j}=δ_{ij}. If the initial state is, however, given in the MCH representation, then c^{MCH}_{j}=δ_{ij} and the determination of c^{diag}=U^{†}c^{MCH} is carried out after the initial QM calculation.
Geometry Input
The initial geometry must be given in a second file in the input format also used by COLUMBUS. The default name for this file is geom. The geometry filename can be given in the input file with the geomfile keyword. Note that the filename has to be enclosed in single or double quotes. See section 4.2 for more details.
Velocity Input
Using the veloc keyword, the initial velocities can be either set to zero, determined randomly or read from a file. Random determination of the velocities is such that each atom has the same kinetic energy, which must be specified after veloc random in units of eV. Determination of the random velocities is detailed in 8.10. Note that after the initial velocities are generated, the RNG is reseeded (i.e., the sequence of random numbers in the surface hopping procedure is independent of whether random initial velocities are used).
Alternatively, the initial velocities can be read from a file.
The default velocity filename is veloc, but the filename can be specified with the velocfile keyword. Note that the filename has to be enclosed in single or double quotes. The file must contain the cartesian components of the velocity for each atom on a new line, in the same order as in the geomety file. The velocity is interpreted in terms of atomic units (bohr/atu). See section 4.3 for more details.
Laser Input
The keyword laser controls whether a laser field is included in the dynamics (influencing the coefficient propagation and the energies/gradients by means of the Stark effect).
The input of an external laser field uses the file laser. This file is specified in 4.5.
Laser Width
In order to detect laserinduced hops, SHARC compares the instantaneous central laser energy with the energy gap between the old and new states. If the difference between the laser energy and the energy gap is smaller than the laser bandwidth, the hop is classified as laserinduced. Those hops are never frustrated and the kinetic energy is not scaled to preserve total energy (instead, the kinetic energy is preserved).
Simulation Timestep
The keyword stepsize controls the length of a timestep (in fs) for the dynamics. The nuclear motion is integrated using the VelocityVerlet algorithm with this timestep. Surface hopping is performed once per timestep and 13 quantum chemistry calculations are performed per timestep (depending on the selection schedule). Each timestep is divided in nsubsteps substeps for the integration of the electronic equationofmotion. Since integration is performed in the MCH representation, the default of 25 substeps is usually sufficient, even if very small potential couplings are encountered.
Simulation Time
The keyword nsteps controls the total length of the simulation. The total simulation time is nsteps times stepsize. nsteps does not include the initial quantum chemistry calculation. Instead of the number of steps the total simulation time can be given directly (in fs) using the keyword tmax. In this case, nsteps is calculated as tmax divided by stepsize. If both keywords (nsteps and tmax) are present, nsteps is used.
Using the keyword killafter, the dynamics can be terminated before the full simulation time. killafter specifies (in fs) the time the trajectory can move in the lowestenergy state before the simulation is terminated. By default, simulations always run to the full simulation time and are not terminated prematurely.
Surface Treatment
The keyword surf controls whether the dynamics runs on diagonal potential energy surfaces (which makes it a SHARC simulation) or on the MCH PESs (which corresponds to a spindiabatic dynamics [24] or FISH [23] simulation). Internally, dynamics on the MCH potentials is conducted by setting the U matrix equal to the unity matrix at each timestep.
Description of Nonadiabatic Coupling
The code allows to propagate the electronic wavefunction using three different quantities describing nonadiabatic effects, see 8.19. The keyword coupling controls which of these quantities is requested from the QM interfaces and used in the propagation. The default is ddr, which requires the nonadiabatic coupling vectors 〈ψ_{α}∂/∂Rψ_{β}〉. For the wavefunction propagation, the scalar product of these vectors and the nuclear velocity is calculated to obtain the matrix 〈ψ_{α}∂/∂tψ_{β}〉. During the propagation, this matrix is interpolated linearly within each classical timestep.
Alternatively, some interfaces can directly calculate the matrix elements 〈ψ_{α}∂/∂tψ_{β}〉, which can be used for the propagation. The corresponding argument to coupling is ddt. In this case, the matrix is taken as constant throughout each classical timestep.
The third possibility is the use of the overlap matrix, requested with coupling overlaps. The overlap matrix is used subsequently in the Local Diabatization algorithm for the wavefunction propagation.
Correction to the Gradients
The correct transformation of the gradients to the diagonal representation includes contributions from the nonadiabatic coupling vectors. Using gradcorrect true, these contributions are included. In this case SHARC will request the calculation of the nonadiabatic coupling vectors, even if they are not used in the wavefunction propagation.
Frustrated Hops and Adjustment of the Kinetic Energy
The keyword ekincorrect controls how the kinetic energy is adjusted after a surface hop to preserve total energy. ekincorrect none deactivates the adjustment, so that the total energy is not preserved after a hop. Using this option, jumps can never be frustrated and are always performed according to the hopping probabilities.
Using ekincorrect parallel_vel, the kinetic energy is adjusted by simply rescaling the nuclear velocities so that the new kinetic energy is E_{tot}−E_{pot}. Jumps are frustrated if the new potential energy would exceed the total energy.
Finally, using ekincorrect parallel_nac, the kinetic energy is adjusted by rescaling the component of the nuclear velocities parallel to the nonadiabatic coupling vector between the old and new state. The hop is frustrated if there is not enough kinetic energy in this direction to conserve total energy. Note that ekincorrect parallel_nac implies the calculation of the nonadiabatic coupling vector, even if they are not used for the wavefunction propagation.
Selection of Gradients and NonAdiabatic Couplings
SHARC allows to selectively calculate only certain gradients and nonadiabatic coupling vectors at each timestep. Those gradients and nonadiabatic coupling vectors not selected are not requested from the interfaces, thus decreasing the computational cost. The selection procedure is detailed in 8.15.
Selection of gradients is activated by grad_select, selection for nonadiabatic couplings by nac_select. Selection is turned off by default.
The selection procedure picks only states which are closer in energy to the classically occupied state than a given threshold. The threshold is 0.5 eV by default and can be adjusted using the eselect keyword.
By default, if SHARC performs such selection it will do two quantum chemistry calls per timestep. In the first call, all quantities are requested except for the ones to be selected. The energies are used to determine which gradients and NACs to calculate in a second quantum chemistry call. The keyword select_directly tells SHARC instead to use the energies of the last timestep, so that only one call per timestep is necessary.
Reference Energy
The keyword ezero gives the energy shift for the diagonal elements of the Hamiltonian. The shift should be chosen so that the shifted diagonal elements are reasonably small (large diagonal elements in the Hamiltonian lead to rapidly changing coefficients, requiring extremely short subtimesteps).
Note that the energy shift default is 0, i.e., SHARC does not choose an energy shift based on the energies at the first time step (this would lead to each trajectory having a different energy shift).
Scaling and Damping
The scaling factor for the energies and gradients must be positive (not zero), see section 8.13.
The damping factor must be in the interval [0,1] (first, since the kinetic energy is always positive; second, because a damping factor larger than 1 would lead to exponentially growing kinetic energy). Also see section 8.3.
RNG Seed
The RNG seed is used to initialize the random number generator, which provides the random numbers for the surface hopping procedure. For details how the seed is used internally, see section 8.14.
Decoherence
The decoherence correction according to Granucci et al. [50] can be applied to the diagonal coefficients by using the keyword decoherence. By default, no decoherence correction is applied. However, the keyword nodecoherence can be used to explicitly turn decoherence off.
The keyword decoherence_param can be used to change the parameter α (see 8.4). The default is 0.1 Hartree, which is the value recommended by Granucci et al. [50].
Ionization
The keyword ionization activates the onthefly calculation of ionization transition properties. If the keyword is given, by default these properties are calculated every timestep. The keyword ionization_step can be used to calculate these properties only every nth timestep.
If the keyword is given, SHARC will request the calculation of the ionization properties from the interface, which needs to be able to calculate them (currently only the COLUMBUS interface allows to perform these calculations in combination with a program which computes the Dyson norms from COLUMBUS output).
Phase Tracking
Using the keywords track_phase and no_track_phase, the tracking of the eigenvector phases of the transformation matrix U
can be switched on or off. However, it is usually not a good idea to deactivate the phase tracking, since this might lead to spurious behaviour in the wavefunction propagation and the surface hopping.
Dipole Moment Gradients
The derivatives of the dipole moments can be included in the gradients. This can be activated with the keyword dipole_gradient. Currently, only the analytical interface can deliver these quantities.
4.1.4 Example
The following input sample shows a typical input for excitedstate dynamics including IC within a singlet manifold plus intersystem crossing to triplet states. It includes a large number of excited singlet states in order to calculate transient absorption spectra. Only the lowest three singlet states actually participate in the dynamics.
nstates 8 0 3 # many singlet states for transient absorption actstates 3 0 3 # only few states to reduce cost stepsize 0.5 # typical timestep for a molecule containing H tmax 1000.0 # one ps surf sharc state 3 mch # start on the S2 singlet state coeff auto # coefficient of S2 will be set to one coupling ddr # \ decoherence #  typical settings ekincorrect parallel_vel # / grad_select # \ nac_select #  improve performance eselect 0.3 # / veloc external # velocity comes from file "veloc" velocfile "veloc" # RNGseed 65435 ezero 399.41494751 # ground state energy of molecule
4.2 Geometry file
The geometry file (default file name is geom) contains the initial coordinates of all atoms. This file must be present when starting a new trajectory.
It uses the COLUMBUS geometry file format (however, ). For each atom, the file contains one line, giving the chemical symbol (a string), the atomic number (a real number), the x, y and z coordinates of the atom in Bohrs (three real numbers), and the relative atomic weight of the atom (a real number). The six items must be separated by spaces. The real numbers are read in using Fortran listdirected I/O, and hence are free format (can have any numbers of decimals, exponential notation, etc.). Element symbols can have at most 2 characters.
The following is an example of a geom file for CH2:
C 6.0 0.0 0.0 0.0 12.000 H 1.0 1.7 0.0 1.2 1.008 H 1.0 1.7 0.0 3.7 1.008
4.3 Velocity file
The velocity file (default veloc) contains the initial nuclear velocities (e.g., from a Wigner distribution sampling). This file is optional (the velocities can be initialized with the veloc input keyword).
The file contains one line of input for each atom, where the order of atoms must be the same as in the geom file. Each line consists of three items, separated by spaces, where the first is the x component of the nuclear velocity, followed by the y and z components (three real numbers). The input is interpreted in atomic units (Bohr/atu).
The following is an example of a veloc file:
0.0001 0.0000 0.0002 0.0002 0.0000 0.0012 0.0003 0.0000 0.0007
4.4 Coefficient file
The coefficient file contains the initial wavefunction coefficients. The file contains one line per state (total number of states, i.e., multiplets count multiple times). Each line specifies the initial coefficient of one state. If the initial state is specified in the MCH representation (input keyword state), then the order of the initial coefficients must be as given by the canonical ordering (see section 8.16). If the initial state is given in diagonal representation, then the initial coefficients correspond to the states given in energetic ordering, starting with the lowest state.
Each line contains two real numbers, giving first the real and then the imaginary part of the initial coefficient of the respective state.
Example:
0.0 0.0 1.0 0.0 0.0 0.0
4.5 Laser file
The laser file contains a table with the amplitude of the laser field ϵ(t) at each timestep of the electronic propagation. Given a laser field of the general form:
each line consists of 8 elements: t (in fs), ℜ(ϵ_{x}(t)), ℑ(ϵ_{x}(t)), ℜ(ϵ_{y}(t)), ℑ(ϵ_{y}(t)), ℜ(ϵ_{z}(t)), ℑ(ϵ_{z}(t)), (all in atomic units), and finally the instantaneous central frequency (also atomic units).
The timestep in the laser file must exactly match the timestep used for the electronic propagation, which is the timestep used for the nuclear propagation (keyword stepsize) divided by the number of substeps (keyword nsubsteps). The first line of the laser file must correspond to t=0 fs.
0.00E+00 0.68E03 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.31E+00 0.10E02 0.77E02 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.33E+00 0.20E02 0.13E01 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.00E+00 0.35E+00 ... ... ... ... ... ... ... ...
Chapter 5 Output files
This chapter documents the content of the output files of SHARC. Those output files are output.log, output.lis, output.dat and output.xyz.
5.1 Log file
The log file output.log contains general information about all steps of the SHARC simulation, e.g., about the parsing of the input files, results of quantum chemistry calls, internal matrices and vectors, etc. The content of the log file can be controlled with the keyword printlevel in the SHARC main input file.
In the following, all printlevels are explained.
Printlevel 0
At printlevel 0, only the execution infos (date, host and working directory at execution start) and build infos (compiler, compile date, building host and working directory) are given.
Printlevel 1
At printlevel 1, also the content of the input file (cleaned of comments and blank lines) is echoed in the log file. Also, the start of each timestep is given.
Printlevel 2
At printlevel 2, the log file also contains information about the parsing of the input files (echoing all enabled options, initial geometry, velocity and coefficients, etc.) and about the initialization of the coefficients after the first quantum chemistry calculation. This printlevel is recommended, since it is the highest printlevel where no output per timestep is written to the log file.
Printlevel 3
This and higher printlevels add output per timestep to the log file. At printlevel 3, the log file contains at each timestep the data from the velocityVerlet algorithm (old and new acceleration, velocity and geometry), the old and new coefficients, the surface hopping probabilities and random number, the occupancies before and after decoherence correction as well as the kinetic, potential and total energies.
Printlevel 4
At printlevel 4, additionally the log file contains information on the quantum chemistry calls (file names, which quantities were read, gradient and nonadiabatic coupling vector selection) and the propagator matrix.
Printlevel 5
At printlevel 5, additionally the log file contains the results of each quantum chemistry calls (all matrices and vectors), all matrices involved in the propagation as well as the matrices involved in the gradient transformation. This is the highest printlevel currently implemented.
5.2 Listing file
The listing file output.lis is a tabular summary of the progress of the dynamics simulation. At the end of each timestep (including the initial timestep), one line with 11 elements is printed. These are, from left to right:
 current step (counting starts at zero for the initial step),
 current simulation time (fs),
 current state in the diagonal representation,
 approximate corresponding MCH state (see subsection 8.11.1),
 kinetic energy (eV),
 potential energy (eV),
 total energy (eV),
 current gradient norm (in eV/Å),
 current expectation values of the state dipole moment (Debye),
 current expectation values of total spin,
 wallclock time needed for the timestep.
The listing file also contains one extra line for each surface hopping event. For accepted hops, the old and new states (in diagonal representation) and the random number are given. Frustrated hops and resonant hops are also mentioned. Note that the extra line for surface hopping occurs before the regular line for the timestep.
The listing file can be plotted with standard tools like GNUPLOT.
Energies
The kinetic energy is calculated at the end of each time step (i.e., after surface hopping events and the corresponding adjustments). The potential energy is the energy of the currently active diagonal state. The total energy is the sum of those two.
Expectation values
The gradient norms given in the listing file is calculated as follows:
which is then transformed to eV/Å.
The expectation values of the dipole moment for the active state β is calculated from:
The expectation value of the total spin of the active state β is calculated as follows:
where S_{α} is the total spin of the MCH state with index α.
5.3 Data file
The data file output.dat contains all relevant data from the simulation for all timesteps, in ASCII format. Accordingly, this file can become quite large for long trajectories or if many states are included, but for most file systems it is easier to deal with a single large file than with many small files.
Usually, after the simulation is finished the data file is processed by data_extractor.x to obtain a number of tabular files which can be plotted. For this, see sections 7.13 for the data extractor and 7.9 for plotting.
5.3.1 Specification of the data file
The data file contains a short header followed by the data per timestep. All quantities are commented in the data file.
The header contains:
 maximum multiplicity,
 number of states per multiplicity,
 number of atoms,
 timestep,
 energy shift,
 flag (overlap matrices),
 flag (laser field),
 number of steps,
 number of substeps,
 full laser field for all substeps (only if flag is set).
The entry for each timestep contains:
 step index
 Hamiltonian in MCH representation,
 transformation matrix U,
 MCH dipole moment matrices (x, y, z),
 overlap matrix in MCH representation (only if flag is set),
 coefficients in the diagonal representation,
 hopping probablities in the diagonal representation
 kinetic energy,
 currently active state in diagonal representation and approximate state in MCH representation,
 random number for surface hopping,
 wallclock time (in seconds)
 geometry (Bohrs),
 velocities (atomic units),
 property matrix.
5.4 XYZ file
The file output.xyz contains the geometries of all timesteps in standard xyz file format. It can be used with visualization programs like MOLDEN, GABEDIT or MOLEKEL to create movies of the molecular motion, or with geo.py (see 7.12) to calculate internal coordinates for each timestep.
The comments of the geometries (given in the second line of each geometry block) contain information about the simulation time and the active state (first in diagonal basis, then in MCH basis).
Chapter 6 Interfaces
This chapter describes the interface between SHARC and quantum chemistry programs. In the first section, the interface is specified (e.g., for users who attempt to create their own interfaces). The description of the currently existing interfaces takes the remainder of this chapter.
6.1 Interface Specifications
From the SHARC point of view, quantum chemical calculation proceeds as follows in the QM directory:
 write a file called QM/QM.in
 call a script called QM/runQM.sh
 read the output from a file called QM/QM.out
For specifications of the formats of these two files (QM.in and QM.out) see below. The executable script QM/runQM.sh must accomplish that all necessary quantum chemical output is available in QM/QM.out.
6.1.1 QM.in Specification
The QM.in file is written by SHARC every time a quantum chemistry calculation is necessary. It contains all information available to SHARC. This information includes the current geometry (and velocity), the timestep, the number of states, the timestep and the unit used to specify the atomic coordinates. The file also contains control keywords and request keywords.
The file format is consistent with a standard xyz file. The first line contains the number of atoms, the second line is a comment. SHARC writes the trajectory ID (a hash of all SHARC input files) to this line. The following lines specify the atom positions. As a fourth, fifth and sixth column, these lines may contain the atomic velocities.
All following lines contain keywords, one per line and possibly with arguments. Comments can be inserted with ‘#’, and empty lines are permissible. Comments and empty lines are only permissible below the xyz file part.
An examplary QM.in file is given in the following:
3 Jobname S 0.0 0.0 0.0 0.000 0.020 0.002 H 0.0 0.9 1.2 0.000 0.030 0.000 H 0.0 0.9 1.2 0.000 0.010 0.000 # This is a comment Init States 3 0 2 Unit Angstrom SOC DM GRAD 1 2 OVERLAP NACDR select 1 2 end
There exist two types of keywords, control keywords and request keywords. Control keywords pass some information to the interface. Request keywords tell the interface to provide a quantity in the QM.out file. Table 6.1 contains all control keywords while table 6.2 lists all request keywords.
Keyword  Description 
init  Specifies that this is the first calculation. The interface should create a save directory to save all information necessary for a restart. 
samestep  Specifies that this is an additional calculation at the same geometry/timestep. 
cleanup  Specifies that all output files of the interface (except QM.out) should be deleted (including the save directory). 
unit  Specifies in which unit the atomic coordinates are to be interpreted. Possible arguments are “angstrom” and “bohr”. 
states  Gives the number of excited states per multiplicity (singlets, doublet, triplets, …). 
dt  Gives the time between the last calculation and the current calculation in atomic units. 
savedir  Gives a path to the directory where the interface should save files needed for restart and between timesteps. If the interfacespecific input files also have this keyword, SHARC assumes that the path in QM.in takes precedence. 
Keyword  Description 
H  Calculate the molecular Hamiltonian (diagonal matrix with the energies of the states of the model space). 
SOC  Calculate the molecular Hamiltonian including the SOCs (not diagonal anymore within the model space). 
DM  Calculate the state dipole moments and transition dipole moments between all states. 
ION  Calculate transition properties between neutral and ionic wavefunctions. 
GRAD  Calculate gradients for all states. If followed by a list of states, calculate only gradients for the specified states. 
NACDT  Calculate the timederivatives 〈Ψ_{1}∂/∂tΨ_{2}〉 by finite differences between the last timestep and the current timestep 
NACDR  Calculate nonadiabatic coupling vectors 〈Ψ_{1}∂/∂RΨ_{2}〉 between all pairs of states. If followed by “select”, read the list of pairs on the following lines until “end” and calculate nonadiabatic coupling vectors between the specified pairs of states. 
OVERLAP  Calculate overlaps 〈Ψ_{1}(t_{0})Ψ_{2}(t)〉 between all pairs of states (between the last and current timestep). If followed by “select”, read the list of pairs on the following lines until “end” and calculate overlaps between the specified pairs of states. 
DMDR  Calculate the cartesian gradients of the dipole moments and transition dipole moments of all states. 
6.1.2 QM.out Specification
The QM.out file communicates back the results of the quantum chemistry calculation to the dynamics code. After SHARC called QM/runQM.sh, it expects that the file QM/QM.out exists and contains the relevant data.
The following quantities are expected in the file (depending whether the corresponding keyword is in the QM.in file): Hamiltonian matrix, dipole matrices, gradients, nonadiabatic couplings (either NACDR or NACDT), overlaps, wavefunction phases, property matrices. The format of QM.out is described in the following.
Each quantity is given as a data block, which has a fixed format. The order of the blocks is arbitrary, and between blocks arbitrary lines can be written. However, within a block no extraneous lines are allowed. Each data block starts with a exclamation mark !, followed by whitespace and an integer flag which specifies the type of data:
1  Hamiltonian matrix 
2  Dipole matrices 
3  Gradients 
4  Nonadiabatic couplings (NACDT) 
5  Nonadiabatic couplings (NACDR) 
6  Overlap matrix 
7  Wavefunction phases(optional) 
8  Wallclock time for QM calculation (optional) 
11  Property matrix (e.g. ionization probabilities) 
12  Dipole moment gradients 
On the next line, two integers are expected giving the dimensions of the following matrix. Note, that all these matrices must be square matrices. On the following lines, the matrix or vector follows. Matrices are in general complex, and real and imaginary part of each element is given as a pair of floating point numbers.
The following shows an example of a 4×4 Hamiltonian matrix. Note that the imaginary parts directly follow the real parts (in this example, the Hamiltonian is real).
! 1 4 4 548.6488 0.0000 0.0000 0.0000 0.0003 0.0000 0.0003 0.0000 0.0000 0.0000 548.6170 0.0000 0.0003 0.0000 0.0003 0.0000 0.0003 0.0000 0.0003 0.0000 548.5986 0.0000 0.0000 0.0000 0.0003 0.0000 0.0003 0.0000 0.0000 0.0000 548.5912 0.0000
The three dipole moment matrices (x, y and z polarization) must follow directly after each other, where the dimension specifier must be present for each matrix. The dipole matrices are also expected to be complexvalued.
! 2 2 2 0.1320 0.0000 0.0020 0.0000 0.0020 0.0000 1.1412 0.0000 2 2 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 2 2 2.1828 0.0000 0.0000 0.0000 0.0000 0.0000 0.6422 0.0000
Gradient and nonadiabatic couplings vectors are written as 3×n_{atom} matrices, with the x, y and z components of one atom per line. These vectors are expected to be real valued. Each vector is preceeded by its dimensions.
6 3 0.0000 6.5429 8.1187 0.0000 5.8586 8.0160 0.0000 6.8428 1.0265 0.0000 6.5429 8.1187 0.0000 5.8586 8.0160 0.0000 6.8428 1.0265
If gradients are requested, SHARC expects every gradient to be present, even if only some gradients are requested. The gradients are expected in the canonical ordering (see section 8.16), which implies that for higher multiplets the same gradient has to be present several times. For example, with 3 singlets and 3 triplets, SHARC expects 12 gradients in the QM.out file (each triplet has three components with M_{s}= 1, 0 or 1).
Similarly, for nonadiabatic coupling vectors, SHARC expects all pairs, even between states of different multiplicity. The vectors are also in canonical ordering, where the inner loop goes over the ket states. For example, with 3 singlets and 3 triplets (12 states), SHARC expects 144 (12^{2}) nonadiabatic coupling vectors in the QM.out file.
! 5 Nonadiabatic couplings (ddr) (2x2x1x3, real) 1 3 ! m1 1 s1 1 ms1 0 m2 1 s2 1 ms2 0 0.0e+0 0.0e+0 0.0e+0 1 3 ! m1 1 s1 1 ms1 0 m2 1 s2 2 ms2 0 +2.0e+0 0.0e+0 0.0e+0 1 3 ! m1 1 s1 2 ms1 0 m2 1 s2 1 ms2 0 2.0e+0 0.0e+0 0.0e+0 1 3 ! m1 1 s1 2 ms1 0 m2 1 s2 2 ms2 0 0.0e+0 0.0e+0 0.0e+0
The nonadiabatic coupling matrix (NACDT keyword), the overlap matrix and the property matrix are single n×n matrices (n is the total number of states), respectively, like the Hamiltonian.
The wavefunction phases are a vector of complex numbers.
The wallclock time is a single real number.
The dipole moment gradients are a list of 3×n_{atom} vectors, each specifying the gradient of one polarization of one dipole moment matrix element. In the outmost loop, the bra index is counted, then the ket index, then the polarization. Hence, the respective entry in QM.out would look like (for 2 states and 1 atom):
! 12 Dipole moment derivatives (2x2x3x1x3, real) 1 3 ! m1 1 s1 1 ms1 0 m2 1 s2 1 ms2 0 pol 0 0.000000000000E+000 0.000000000000E+000 0.000000000000E+000 1 3 ! m1 1 s1 1 ms1 0 m2 1 s2 1 ms2 0 pol 1 0.000000000000E+000 0.000000000000E+000 0.000000000000E+000 1 3 ! m1 1 s1 1 ms1 0 m2 1 s2 1 ms2 0 pol 2 0.000000000000E+000 0.000000000000E+000 0.000000000000E+000 1 3 ! m1 1 s1 1 ms1 0 m2 1 s2 2 ms2 0 pol 0 1.000000000000E+000 0.000000000000E+000 0.000000000000E+000 1 3 ! m1 1 s1 1 ms1 0 m2 1 s2 2 ms2 0 pol 1 0.000000000000E+000 0.000000000000E+000 0.000000000000E+000 1 3 ! m1 1 s1 1 ms1 0 m2 1 s2 2 ms2 0 pol 2 0.000000000000E+000 0.000000000000E+000 0.000000000000E+000 1 3 ! m1 1 s1 2 ms1 0 m2 1 s2 1 ms2 0 pol 0 1.000000000000E+000 0.000000000000E+000 0.000000000000E+000 1 3 ! m1 1 s1 2 ms1 0 m2 1 s2 1 ms2 0 pol 1 0.000000000000E+000 0.000000000000E+000 0.000000000000E+000 1 3 ! m1 1 s1 2 ms1 0 m2 1 s2 1 ms2 0 pol 2 0.000000000000E+000 0.000000000000E+000 0.000000000000E+000 1 3 ! m1 1 s1 2 ms1 0 m2 1 s2 2 ms2 0 pol 0 0.000000000000E+000 0.000000000000E+000 0.000000000000E+000 1 3 ! m1 1 s1 2 ms1 0 m2 1 s2 2 ms2 0 pol 1 0.000000000000E+000 0.000000000000E+000 0.000000000000E+000 1 3 ! m1 1 s1 2 ms1 0 m2 1 s2 2 ms2 0 pol 2 0.000000000000E+000 0.000000000000E+000 0.000000000000E+000
6.1.3 Further Specifications
The interfaces may require additional input files beyond QM.in, which contain static information. This may include paths to the executable quantum chemistry program, paths to scratch directories or input templates for the quantum chemistry calculation (e.g. active space specifications, basis sets, etc.).
The dynamics code does not depend on these additional files.
6.1.4 Save Directory Specification
The interfaces must be able to save all information necessary for restart to a given directory. The absolute path is written to QM.in by SHARC. Hence, for the trajectories the path to the save directory is always a subdirectory of the working directory of SHARC.
6.2 MOLPRO Interface
The SHARC– MOLPRO interface allows to run SHARC dynamics with MOLPRO‘s CASSCF wavefunctions. RASSCF is not supported, since on RASSCF level stateaveraging over different multiplicities is not possible. The interface uses MOLPRO‘s CI program in order to calculate transition dipole moments and spinorbit couplings. Gradients and nonadiabatic coupling vectors are calculated using MOLPRO‘s CADPACK code (hence no generally contracted basis sets can be used). Time derivatives and overlaps are obtained with the DDR prodecure. Hence, all types of couplings usable by SHARC are available through this interface. The interface also allows to calculate the expectation value of the angular momentum operators.
The SHARC– MOLPRO interface needs two additional input files, which should be present in QM/. Those input files are SH2PRO.inp, which contains the paths to MOLPRO and the scratch directory, and MOLPRO.template, which is a minimal input from which the full input can be built. If QM/wf.init is present, it will be used as a MOLPRO wavefunction file containing the initial MOs.
6.2.1 Interface specific input file: SH2PRO.inp
The interface requires some additional information beyond the content of QM.in. This information is given in the file SH2PRO.inp, which must reside in the directory where the interface is started. This file uses a simple “keyword argument” syntax. Comments using # and blank lines are possible, the order of keywords is arbitrary. Lines with unknown keywords are ignored, since the interface just searches the file for certain keywords.
Table 6.3 lists the existing keywords.
Keyword  Description 
molpro  Is followed by a string giving the path to the MOLPRO executable. Relative and absolute paths, environment variables and can be used. 
scratchdir  Is a path to the temporary directory. Relative and absolute paths, environment variables and can be used. If the directory does not exist, the interface will create it. In any case, the interface will delete this directory after the calculation. 
gradaccudefault  (float) Default accuracy for CPMCSCF. 
gradaccumax  (float) Worst acceptable accuracy for CPMCSCF (see below). 
checknacs  (boolean) Use the MRCI overlaps to check whether the time derivatives and overlap matrices are correct. 
correctnacs  (boolean) Replace bad time derivatives with the scalar product of nonadiabatic coupling vectors and velocities. 
checknacs_mrcio  (float) Threshold for the MRCI overlaps to be considered bad. 
checknacs_ediff  (float, eV) If MRCI overlaps indicate inaccurate couplings, set to zero couplings between states which are farther apart than this value. 
Note that only the first two keywords are mandatory. The last four keywords are only used for calculations with time derivatives or overlaps. The checking and correction of these couplings has to be considered experimental. Normally, checknacs should be “False” (in this case the other three keywords are not necessary).
6.2.2 Template file: MOLPRO.template
The template file is a MOLPRO input file specifying a stateaveraged CASSCF calculation. There must not be any file specifications and no geometry input in this file. It should contain memory specifications, basis set, and optionally settings like DouglasKroll.
Most importantly, it has to contain a CASSCF input block with the number of frozen, closedshell and occupied orbitals. Additionally, all wavefunctions for the stateaveraging have to be defined. For each multiplicity, the wf, state and weight keywords must be present. The following shows an example input for CASSCF(12,9) with 4 singlets and 3 triplets in the stateaveraging.
***,Example memory,100,M dkroll=1 dkho=2 basis=631G* {casscf frozen,0 closed,23 occ,32 wf,58,1,0 state,4 weight,1,1,1,1 wf,58,1,2 state,3 weight,1,1,1 };
6.2.3 Error checking
The interface is written such that the output of MOLPRO is checked for commonly occuring errors, mostly bad convergence in the MCSCF or CPMCSCF parts. In these cases, the input is adjusted and MOLPRO restarted. This will be done until all calculations are finished or an unrecoverable error is detected.
The interface will try to solve the following error messages:
EXCESSIVE GRADIENT IN CI This error message can occur in the MCSCF part. The calculation is restarted with a Pspace threshold (see MOLPRO manual) of 1. If the error remains, the threshold is quadrupled until the calculation converges or the threshold is above 100.
NO CONVERGENCE IN REFERENCE CI The error occurs in the CI part. The calculation is restarted with a Pspace threshold (see MOLPRO manual) of 1. If the error remains, the threshold is quadrupled until the calculation converges or the threshold is above 100.
NO CONVERGENCE OF CPMCSCF This error occurs when solving the linear equations needed for the calculation of MCSCF gradients or nonadiabatic coupling vectors. In this case, the interface finds in the output the value of closest convergence and restarts the calculation with the value found as the new convergence criterion. This ensures that the CPMCSCF calculation converges, albeit with lower accuracy for this gradient for this timestep.
This error check is controlled by two keywords in the SH2PRO.inp file. The interface first tries to converge the CPMCSCF calculation to gradaccudefault. If this fails, it tries to converge to the best value possible within 900 iterations. gradaccumax defines the worst accuracy accepted by the interface. If a CPMCSCF calculation cannot be converged below gradaccumax then the interface exits with an error, leading to the abortion of the trajectory.
6.2.4 Things to keep in mind
Initial orbital guess
For CASSCF calculations it is always a good idea to start from converged MOs from a nearby geometry. For the first timestep, if a file QM/wf.init is present, the SHARC– MOLPRO interface will take this file for the starting orbitals. In subsequent calculations, the MOs from the previous step will be used.
CPMCSCF calculations vs. DDR calculations
In the current version of the interface, the keywords grad and nacdr cannot be combined with the keywords nacdt and overlap. This is related to the way the interface allocates the records for gradients in the wavefunction file of Molpro. This issue might be resolved in future releases.
This has the consequence that for trajectories using time derivatives or overlaps the gradient selection must be activated (though the selection criterion can be chosen arbitrarily high). Also, the select_directly keyword is not compatible with time derivatives and overlaps. setup_traj.py automatically takes care of this constraints.
Basis sets
Note that Molpro cannot calculation SACASSCF gradients for generally contracted basis sets (like Dunning’s “cc” basis sets or Roos’ “ANO” basis sets). Only segmented basis sets are allowed, like the Pople basis sets and the “def” family from TURBOMOLE.
6.2.5 Molpro input generator: molpro_input.py
In order to quickly setup simple calculations using MOLPRO, the SHARC suite contains a small script called molpro_input.py. It can be used to setup singlepoint calculations, optimizations and frequency calculations on the HF, DFT, MP2 and CASSCF level of theory. Of course, MOLPRO has far more capabilities, but these are not covered by molpro_input.py. However, molpro_input.py can also prepare template files which are compatible with the SHARC–Molpro interface.
The script interactively asks the user to specify the calculation and afterwards writes an input file and optionally a run script.
Input
Type of calculation
Choose to either perform a singlepoint calculation or an optimization (including optionally frequency calculation), or to generate a template file. In the latter case, no geometry file is needed. The script looks for a MOLPRO.input in the same directory and allows to copy the settings.
For singlepoint calculations, optimizations and frequency calculations, files in MOLDEN format called geom.molden, opt.molden or freq.molden, respectively, are created (containing the orbitals, optimization steps and normal modes, respectively). The file freq.molden can be used to generate initial conditions with wigner.py.
Geometry
Specify the geometry file in xyz format. Number of atoms and total nuclear charge is detected automatically. After the user inputs the total charge, the number of electrons is calculated automatically.
In the case of the generation of a template file, instead only the number of electrons is required.
Nondefault atomic masses
If a frequency calculation is requested, the user may modify the mass of specific atoms (e.g. to investigate isotopic effects). In the following menu, the user can add or remove atoms with their mass to a list containing all atoms with nondefault masses. Each atom is referred to by its number as in the geometry file. Using the command show the user can display the list of atoms with nondefault masses. Typing end confirms the list.
Note that when using the produced MOLDEN file later with wigner.py, the user has to enter the same nondefault masses again, since the MOLDEN file does not contain the masses and wigner.py has no way to retrieve these numbers.
Level of theory
Supported are HartreeFock (HF), Density Functional Theory (DFT), MøllerPlesset perturbation theory (MP2) and CASSCF (either singlestate or stateaveraged). All methods are compatible with oddelectron wavefunctions (molpro_input.py will use the corresponding UHF, UMP2 and UKS keywords in the input file, if necessary).
For template generation stateaverage CASSCF is automatically chosen. All methods can be combined with optimizations and frequency calculations, however, the frequency calculation is much more efficient with HF or SSCASSCF.
DFT functional
For DFT calculations, enter a functional and choose whether dispersion correction should be applied. Note that the functional is just a string which is not checked by molpro_input.py.
Basis set
The basis set is just a string which is not checked by molpro_input.py.
CASSCF settings
For CASSCF calculations, enter the number of active electrons and orbitals.
For SSCASSCF, only the multiplicity needs to be specified. For SACASSCF, specify the number of states per multiplicity to be included. Note that MOLPRO allows to average over states with different numbers of electrons. This feature is not supported in molpro_input.py. However, the user can generate a closelymatching input and simply add the missing states to the CASSCF block manually.
For optimizations at SACASSCF level, the state to be optimized has to be given.
Memory
Enter the amount of memory for MOLPRO. Note that values smaller than 50 MB are ignored, and 50 MB are used in this case.
Run script
If requested, the script also generates a simple Bash script (run_molpro.sh) to directly execute MOLPRO. The user has to enter the path to MOLPRO and the path to a suitable (fast) scratch directory.
Note that the scratch directory will be deleted after the calculation, only the wavefunction file wf will be copied back to the main directory.
6.3 MOLCAS Interface
The SHARC– MOLCAS interface can be used to conduct excitedstate dynamics based on MOLCAS‘ CASSCF wavefunctions. RASSCF wavefunctions are not supported currently. The interface uses the modules GATEWAY, SEWARD (integrals), RASSCF (wavefunction, energies), RASSI (transition dipole moments, spinorbit couplings, overlaps), MCLR and ALASKA (gradients). Since MOLCAS is not able to calculate the full nonadiabatic coupling vectors, only wavefunction overlaps are currently implemented in the interface.
The interface needs two additional input files, a template file for the quantum chemistry (file name is MOLCAS.template) and a general input file (SH2CAS.inp). If QM/MOLCAS.<i>.JobIph.old files are present, they are used as initial wavefunction files, where <i> is the multiplicity.
6.3.1 Interface specific input file: SH2CAS.inp
The file SH2CAS.inp contains mainly paths (to the MOLCAS executables, to the scratch directory, etc.). This file must reside in the same directory where the interface is started. It uses a simple “keyword argument” syntax. Comments using # and blank lines are possible, the order of keywords is arbitrary. Lines with unknown keywords are ignored, since the interface just searches the file for certain keywords.
Table 6.4 lists the existing keywords.
Keyword  Description 
molcas  Is the path to the MOLCAS installation. Relative and absolute paths, environment variables and can be used. The interface will set $MOLCAS to this path. If this keyword is not present in SH2CAS.inp, the interface will use the environment variable $MOLCAS, if it is set. 
scratchdir  Is a path to the temporary directory. Relative and absolute paths, environment variables and can be used. If it does not exist, the interface will create it. In any case, the interface will delete this directory after the calculation. 
memory  The memory usable by MOLCAS. The interface will set $MOLCASMEM to this value. The default is 10 MB. 
project  Can be used to set a MOLCAS project name ($Project). If not used, the interface will generate a default project name. 
Note that the interface sets all environment variables necessary to run MOLCAS (e.g., $MOLCAS, $MOLCASMEM, $WorkDir ,$Project) automatically, based on the paths to MOLCAS and the scratch directory from SH2CAS.inp.
6.3.2 Template file: MOLCAS.template
This file contains the specifications for the wavefunction. Note that this is not a valid MOLCAS input file. No sections like $GATEWAY, etc., can be used. The file only contains a number of keywords, given in table 6.5.
Keyword  Description 
basis  The basis set used. Note that some basis sets (e.g., Pople basis sets) do not work, since the spinorbit integrals cannot be calculated. 
ras2  Number of active orbitals for CASSCF. 
nactel  Number of active electrons for CASSCF. 
inactive  Number of inactive orbitals. 
spin  The full line reads as spin s roots r and specifies the number of roots for the given multiplicity. This influences only the number of states in the stateaveraging procedure. The number of states in the dynamics must not be larger than the number of states given here. 
The template file can be setup with the tool molcas_input.py.
6.3.3 Template file generator: molcas_input.py
This is a small interactive script to generate template files for the SHARC– MOLCAS interface. It simply queries the user for some input parameters and then writes the file SH2CAS.inp, which can be used to run SHARC simulations with the SHARC– MOLCAS interface.
Geometry file
The geometry file is only used to calculate the nuclear charge.
Charge
This is the overall charge of the molecule. This number is used with the nuclear charge to calculate the number of electrons.
Basis set
This is simply a string, which is not checked by the script to be a valid basis set of the MOLCAS library.
Number of active electrons and orbitals
These settings are necessary for the definition of the CASSCF wavefunction. The number of inactive orbitals is automatically calculated from the total number of electrons and the number of active electrons.
States for stateaveraging
For each multiplicity, the number of states for the stateaveraging procedure must be equal or larger than the number of states used in the dynamics.
6.4 COLUMBUS Interface
The SHARC– COLUMBUS interface allows to run SHARC dynamics based on COLUMBUS‘ CASSCF, RASSCF, MRCI and LRTMRAQCC wavefunctions. The interface is compatible to COLUMBUS calculations utilizing the COLUMBUS– MOLCAS interface only ( SEWARD integrals and ALASKA gradients). Unfortunately, time derivatives and nonadiabatic coupling vectors cannot be calculated with this setup, only wavefunction overlaps can be obtained. The interface utilizes the cioverlap program by Jiri Pittner [27] to calculate the overlap matrices needed for local diabatization propagation. The interface can also calculate Dyson norms between neutral and ionic wavefunctions using a suitable usersupplied code.
The interface needs as additional input the file QM/SH2COL.inp and a template directory containing all input files needed for the COLUMBUS calculations. Additionally, initial MOs can be given in the file QM/mocoef_mc.init.
6.4.1 Template input
The interface does not generate the full COLUMBUS input onthefly. Instead, the interface uses an existing set of input files and performs only necessary modifications (e.g., the number of states). The set of input files must be provided by the user. Please see the COLUMBUS online documentation and, most importantly, the
COLUMBUS SOCI tutorial for a documentation of the necessary input.
Generally, the input consists of a directory with one subdirectory with input for each multiplicity (singlets, doublets, triplets, …). However, evenelectron wavefunctions of different multiplicities can be computed together in the same job if spinorbit couplings are desired. Otherwise independent multipleDRT inputs (ISC keyword) are also acceptable. Note that symmetry is not allowed when using the interface.
The path to the template directory must be given in SH2COL.inp.
Integral input
The interface is only able to use input for calculations using MOLPRO/SEWARD integrals. Also make sure that you include scalarrelativistic (DouglasKrollHess) and spinorbit (AMFI) effects in the integral input.
It is important to make sure that the order of atoms in the template input files and in the SHARC input is consistent.
MCSCF input
The MCSCF section can use any desired stateaveraging scheme. However, frozen core orbitals in the MCSCF step are not possible (since otherwise gradients cannot be computed). Prepare the MCSCF input for CI gradients. It is advisable to use very tight MCSCF convergence criteria.
MRCI input
Either prepare a singleDRT input without SOCI (to cover a single multiplicity), a singleDRT input with SOCI and a sufficient maximum multiplicity for spinorbit couplings or a independent multipleDRT input (ISC case). Make sure that all multiplicities are covered with all input directories.
In the MRCI input, make sure to use sequential ciudg. Also take care to setup gradient input on MRCI level.
Job control
Setup a singlepoint calculation with the following steps:
 SCF
 MCSCF
 MRCISD (serial operation) or SOCI coupled to nonrel CI (for SOCI DRT inputs)
 oneelectron properties for all methods
 transition moments for MRCISD
 nonadiabatic couplings (and/or gradients)
Request first transition moments and interstate couplings in the following dialogues. Analysis in internal coordinates and intersection slope analysis are not required.
6.4.2 Interface specific input file: SH2COL.inp
Beyond the information from QM.in the interface needs additional input, which is read from the file SH2COL.inp. Table 6.6 lists the keywords which can be given in the file.
Keyword  Description 
columbus  Is followed by a string giving the path to the COLUMBUS main directory. Relative and absolute paths, environment variables and can be used. 
molcas  Is followed by a string giving the path to the MOLCAS main directory. Relative and absolute paths, environment variables and can be used. This path is only used for overlap calculations (since in this case the interface calls MOLCAS explicitly). Otherwise COLUMBUS will use the path to MOLCAS specified during the installation of COLUMBUS. 
scratchdir  Is a path to the temporary directory. Relative and absolute paths, environment variables and can be used. If it does not exist, the interface will create it. In any case, the interface will delete this directory after the calculation. 
savedir  Is a path to another temporary directory. Relative and absolute paths, environment variables and can be used. The interface will store files needed for restart there. 
dyson  Path to the dyson code. Relative and absolute paths, environment variables and can be used. Only necessary if Dyson norms are calculated. Needs a suitable code that computes Dyson norms. 
civecconsolidate  Path to the civecconsolidate code suitable for preparation of input for the Dyson program (CI vectors in terms of determinants). Relative and absolute paths, environment variables and can be used. Only used if Dyson norms are calculated. If not present, the civecconsolidate code in the $COLUMBUS directory is used. Note that in the current version of COLUMBUS (Oct. 2014) civecconsolidate is not compatible with the calculation of Dyson norms. 
memory  (integer, MB) Memory for COLUMBUS. Note that the maximum amount of memory used by cioverlaps cannot be controlled with this keyword. 
ciothres  (float) Threshold for including a pair of determinants in the cioverlaps program. 
dysonthres  (float) Threshold for including a pair of determinants in the calculation of Dyson norms. 
nooverlap  (no argument) Do not keep the binary files necessary to calculate wavefunction overlaps in the next timestep. 
template  Is followed by the path to the directory containing the template subdirectories. Relative and absolute paths, environment variables and can be used. 
DIR  See 6.4.3. 
MOCOEF  See 6.4.3. 
6.4.3 Template setup
The template directory contains several subdirectories with input for different multiplicities. An example is given in figure 6.1. In SH2COL.inp, the user has to associate each multiplicity to a subdirectory. The line “DIR 1 Sing_Trip” would make the interface use the input files from the subdirectory Sing_Trip when calculating singlet states (the 1 refers to singlet calculations). All calculations using a particular input subdirectory are called a job.
Additionally, the user must specify which job(s) provide the MO coefficients (e.g., the calculation for doublet states could be based on the same MOs as the singlet and triplet calculation). The line “MOCOEF Doub_Quar Sing_Trip” would tell the interface to do a MCSCF calculation in the Sing_Trip job, and reuse the MOs when doing the Doub_Quar job without reoptimizing the MOs.
6.5 Analytical PESs Interface
The SHARC suite also contains an interface which allows to run dynamics simulations on PESs expressed with analytical functions.
In order to allow dynamics simulations in the same way as it is done onthefly, the interface uses two kinds of potential couplings:
 couplings that are prediagonalized in the interface (yielding the equivalent of the MCH basis),
 couplings given by the interface to SHARC as offdiagonal elements.
The interface needs one additional input file, called SH2Ana.inp, which contains the definitions of all analytical expressions.
6.5.1 Parametrization
The interface has to be provided with analytical expressions for all matrix elements of the following matrices in the diabatic basis:
 Hamiltonian: H
 Derivative of the Hamiltonian with respect to each atomic coordinate: H_{xi}
 (Transition) dipole matrices for each polarization direction: M_{i}
 Real and imaginary part of the SOC matrix: Σ
 (Optionally) the derivatives of the transition dipole matrices: D_{i,xj}
The diabatic Hamiltonian is diagonalized:
Then the following calculations lead to the MCH matrices which are passed to SHARC:
The MCH Hamiltonian is the diagonalized diabatic Hamiltonian plus the SO matrix transformed into the MCH basis. The gradients in the MCH basis are obtained by transforming the derivative matrices into the MCH basis. The dipole matrices are also simply transformed into the MCH basis. The overlap matrix is the overlap of old and new transformation matrix.
6.5.2 Interfacespecific input file
The interfacespecific input file is called SH2Ana.inp. It contains the analytical expressions for all matrix elements mentioned above. All analytical expressions in this file are evaluated considering the atomic coordinates read from QM.in.
The file consists of a file header and the file body. The file body consists of variable definition blocks and matrix blocks.
Header
The header looks similar to an xyz file:
2 2 I xI 0 0 Br xBr 0 0
Here, the first line gives the number of atoms and the second line the number of states.
On the remaining lines, each cartesian component of the atomic coordinates is associated to a variable name, which can be used in the analytical expressions. If a zero (0) is given instead of a variable name, then the corresponding cartesian coordinate is neglected. In the above example, the variable name xI is associated with the x coordinate of the first atom given in QM.in. The y and z coordinates of the first atom are neglected.
All variable names must be valid Python identifiers and must not start with an underscore. Hence, all strings starting with a letter, followed by an arbitrary number of letters, digits and underscores, are valid. It is not allowed to use a variable name twice.
Note that the file header also contains the atom labels, which are just used for crosschecking against the atom labels in QM.in.
The file header must not contain comments, neither at the end of a line nor separate lines. Also, blank lines are not allowed in the header. After the last line of the header (where the variables for the n_{atom}th atom are defined), blank lines and comments can be used freely (except in matrix blocks).
Variable definition blocks
Variable definition blocks can be used to store additional numerical values (beyond the atomic coordinates) in variables, which can then be used in the equations in the matrix blocks. The most obvious use for this is of course to define values which will appear several times in the equations.
A variable definition block looks like:
Variables A1 0.067 g1 0.996 # Trailing comment # Blank line with comment only R1 4.666 End
Each block starts with the keyword “Variables” and is terminated with “End”. Inbetween, on each line a variable name and the corresponding numerical value (separated by blanks) can be given. Note that the naming conventions given above also apply to variables defined in these blocks.
There can be any number of complete variable definitions blocks in the input file. All blocks are read first, before any matrix expressions are evaluated. Hence, the relative order of the variable blocks and the matrix blocks does not matter. Also, note that variable names must not appear twice, so variables cannot be redefined halfway through the file.
Matrix blocks
The most important information in the input file are of course contained in the expressions in the matrix blocks. In general, a matrix block has the following format:
Matrix_Identifier V11 V12, V22 V13, V23, V33 ...
The first line identifies the type of matrix. Those are valid identifiers:
Hamiltonian  Defines the Hamiltonian including the diabatic potential couplings. 
Derivatives followed by a variable name  Derivative of the Hamiltonian with respect to the given variable. 
Dipole followed by 1, 2 or 3  (Transition) dipole moment matrix for cartesian direction x, y or z, respectively. 
Dipolederivatives followed by 1, 2 or 3 followed by a variable name  Derivative of the respective dipole moment matrix. 
SpinOrbit followed by R or I  Real or Imaginary (respectively) part of the spinorbit coupling matrix Σ. 
Since the interface searches the file for these identifiers starting from the top until it is found, for each matrix only the first block takes effect. Note that the Hamiltonian and all relevant derivatives must be present. If dipole matrix, dipole derivative matrix or SO matrix definitions are missing, they will be assumed zero.
In the lines after the identifier, the expressions for each matrix element are given. Note the lower triangular format (all matrices are assumed symmetric, except the imaginary part of the SO matrix, which is assumed antisymmetric). Matrix elements must be separated by commas (so that whitespace can be used inside the expressions). There must be at least as many lines as the number of states (additional lines are neglected). If any line or matrix element is missing, the interface will abort.
An examplary block looks like:
Hamiltonian A1*( (1.math.exp(g1*(R1xI+xBr)))**21.), 0.0006, 3e5*(xIxBr)**2
It is important to understand that the expressions are directly evaluated by the Python interpreter, hence all expressions must be valid Python expressions which evaluate to numeric (integer or float) values. Only the variables defined above can be used.
Note that exponentiation in Python is **. In order to provide most usual mathematical functions, the math module is available. Among others, the math module provides the following functions:
 math.exp(x): Exponential function
 math.log(x): Natural logarithm
 math.pow(x,y): x^{y}
 math.sqrt(x): √x
 math.cos(x), math.sin(x), math.tan(x)
 math.acos(x), math.asin(x), math.atan(x)
 math.atan2(y,x): tan^{−1}([y/x]), as in many programming languages (takes care of phases)
 math.pi, math.e: π and Euler’s number
 math.cosh(x), math.sinh(x): Hyperbolic functions (also tanh, acosh, asinh, atanh are available)
Chapter 7 Auxilliary Scripts
In this chapter, all auxiliary scripts and programs are documented. Input generators (currently molpro_input.py and molcas_input.py) are documented in the relevant interface sections.
All auxiliary scripts are either interactive – prompting user input from stdin in order to setup a certain task – or noninteractive, meaning they are controlled by commandline arguments and options, in the same way as many commandline tools work.
All interactive scripts sequentially ask a number of questions to the user. In many cases, a default value is presented, which is either preset or detected by the scripts based on the availability of certain files. Furthermore, the scripts feature autocompletion of paths and filenames (use TAB), which is active only in questions where autocompletion is relevant.
All interactive scripts write a file called KEYSTROKES.<script_name> which contains the user input from the last completed session. These files can be piped to the interactive scripts to perform the same task again, for example:
user@host> cat KEYSTROKES.excite   $SHARC/excite.py
Note the , which tells cat to switch to stdin after the file ends, so that the user can proceed if the script asks for more input than contained in the KEYSTROKES file.
All noninteractive scripts can be called with the h option to obtain a short description, usage information and a list of the command line options.
All scripts can be safely killed during a run by using CTRL+C. In the case of interactive scripts, a KEYSTROKES.tmp file remains, containing the user input made so far.
7.1 Wigner Distribution Sampling: wigner.py
The first step in preparing the dynamics calculation is to obtain a set of physically reasonable initial conditions. Each initial condition is a set of initial atomic coordinates, initial atomic velocities and initial electronic state. The initial geometry and velocities can be obtained in different ways. With SHARC, usually sampling of a quantumharmonic Wigner distribution is performed.
The sampling is carried out with the noninteractive Python script wigner.py. The theoretical background is summarized in section 8.12.
7.1.1 Usage
The general usage is
user@host> python $SHARC/wigner.py [options] filename.molden
wigner.py takes exactly one commandline argument (the input file with the frequencies and normal modes), plus some options. Usually, the n option is necessary, since by default only 3 initial conditions are created.
The argument is the filename of the file containing the information about the vibrational frequencies and normal modes. The file is by default assumed to be in the MOLDEN format. For usage with wigner.py, only the following blocks have to be present:

Alternatively, the information about the vibrational frequencies and normal modes can be read directly from a MOLPRO output file of a frequency calculation. In this case, in addition to the filename also the switch M has to be given:
user@host> python $SHARC/wigner.py M MOLPRO.out
The script accepts a number of commandline options, specified in table 7.1.
Option  Description  Default 
h  Display help message and quit.  – 
m  Modify atom masses  Most common isotope. 
M  Assume a MOLPRO input file.  Assume a MOLDEN file. 
n INTEGER  Number of initial conditions to generate.  3 
o FILENAME  Output filename.  initconds 
r INTEGER  Seed for random number generator.  16661 
s FLOAT  Scaling factor for the frequencies.  1.0 
7.1.2 Output
The script wigner.py generates a single output file, by default called initconds. All information about the initial conditions is stored in this file. Later steps in the preparation of the initial conditions add information about the excited states to this file. The file is formatted in a humanreadable form.
The initconds file format is specified in section 7.3.5.
7.1.3 Nondefault Masses
When the m option is used, the script will ask the user to interactively modify the atom masses. For each atom (referred to by the atom index as in the MOLDEN file), a mass can be given (relative atomic weights). Note that the frequency calculation which produces the MOLDEN or MOLPRO file should be done with the same atomic masses.
7.2 Setup of Initial Calculations: setup_init.py
The interactive script setup_init.py creates input for singlepoint calculations at the initial geometries given in an initconds file. These calculations might be necessary for some schemes to select the initial electronic state of the trajectory, based on the excitation energies and oscillator strength of the transitions from ground state to the excited state under consideration.
There are other choices of the initial state possible, which do not require singlepoint calculations at all initial geometries. See the description of excite.py (section 7.3). In this case, setup_init.py can be used to set up only the calculation at the equilibrium geometry (see below at “Range of Initial Conditions”).
7.2.1 Usage
The script is interactive, and can be started by simply typing
user@host> python $SHARC/setup_init.py
Please be aware that the script will setup the calculations in the directory where it was started, so the user should cd to the desired directory before executing the script.
Please note that the script does not expand or shell variables, except where noted otherwise.
7.2.2 Input
The script will prompt the user for the input. In the following, all input parameters are documented:
Initial Conditions File
Enter the filename of the initial conditions file, which was generated beforehand with wigner.py. If the script finds a file called initconds, the user is asked whether to use this file, otherwise the user has to enter an appropriate filename. The script detects the number of initial conditions and number of atoms automatically from the initial conditions file.
Range of Initial Conditions
The initial conditions in initconds are indexed, starting with the index 1. In order to prepare ab initio calculations for a subset of all initial conditions, enter a range of indices, e.g. a and b. This will prepare all initial conditions with indices in the interval [a,b]. In any case, the script will additionally prepare a calculation for the equilibrium geometry (except if a finished calculation for the equilibrium geometry was found).
If the interval [0,0] is given, the script will only setup the calculation at the equilibrium geometry.
Number of states
Here the user can specify the number of excited states to be calculated. Note that the ground state has to be counted as well, e.g. if 4 singlet states are specified, the calculation will involve the S_{0}, S_{1}, S_{2} and S_{3}. Also states of higher multiplicity can be given, e.g. triplet or quintett states.
For evenelectron molecules, including oddelectron states (e.g. doublets) is only useful if transition properties for ionization can be computed (e.g. Dyson norms with the COLUMBUS interface). These transition properties can be used to calculate ionization spectra or to obtain initial conditions for dynamics after ionization.
Spinorbit calculation
Usually it is sufficient to calculate the spinfree excitation energies and oscillator strengths in order to decide for the initial state. However, using this option, the effects of spinorbit coupling on the excitation energies and oscillator strengths can be included. Note that the script will never calculate spinorbit couplings if only singlet states are included in the calculation.
Interface
In this point, choose any of the displayed interfaces to carry out the ab initio calculations. Enter the corresponding number.
The following input section depends on the chosen interface.
7.2.3 Input for MOLPRO
Path to MOLPRO
Here the user is prompted to provide the path to the MOLPRO executable.
Note that the setup script will not expand the user ( ) and shell variables (since the calculations might be running on a different machine than the one used for setup, so the meaning of shell variables could be different). and shell variables will only be expanded by the interfaces during the actual calculation.
Scratch directory
The script takes the string without any checking. Each individual initial condition uses a subdirectory ICOND_%05i with the appropriate number. and shell variables will only be expanded by the interface during the actual calculation.
Note that you can use, e.g., ./temp/ as scratch directory, which is a subdirectory of the working directory of the interface. Using ./ as scratch directory is not recommended, since the interface will delete the scratch directory.
Template file
Enter the filename for the MOLPRO template input. This file contains the details of the ab initio calculation, like basis set, active orbitals and electrons, number of electrons and number of states for stateaveraging. For MOLPRO, it also contains the memory usage. For details, see the section about the SHARC– MOLPRO interface (6.2). The setup script will check whether the template file contains the necessary entries.
Initial wavefunction file
You can provide an initial wavefunction (with an appropriate MO guess) to MOLPRO. This usually provides a drastic speedup for the CASSCF calculations and helps in ensuring that all calculations are based on the same active space. In simple cases it might be acceptable to not provide an initial wavefunction.
7.2.4 Input for COLUMBUS
Path to COLUMBUS
Here the user is prompted to provide the path to the COLUMBUS directory. Note that the script will not expand the user ( ) and shell variables (since possibly the ab initio calculations are running on a different machine than the one used for setup). and shell variables will only be expanded by the interfaces during the actual calculation.
Note that $COLUMBUS does not have to be defined on the setup host, but it has to be defined on the hosts running the calculations.
Scratch directory
The script takes the string without any checking. Each individual initial condition uses a subdirectory ICOND_%05i with the appropriate number. and shell variables will only be expanded by the interface during the actual calculation.
Note that you can use, e.g., ./temp/ as scratch directory, which is a subdirectory of the working directory of the interface. Using ./ as scratch directory is not recommended, since the interface will delete the scratch directory.
Template directory
Enter the path to a directory containing subdirectories with COLUMBUS input files necessary for the calculations. The setup script will expand and shell variables and will pass the absolute path to the calculations.
The script will autodetect (based on the cidrtin files) which subdirectory contains the input for each multiplicity. The user has to check in the following dialog whether the association of multiplicities with job directories is correct. Additionally, the association of MO coefficient files to the jobs has to be checked. See the section on the SHARC– COLUMBUS interface (6.4) for further details.
Note that the setup script does not check any content of the input files beyond the multiplicity. Note that transmomin and control.run do not need to be present (and that their content has no effect on the calculation), since these files are written onthefly by the interface.
Initial orbital coefficient file
You can provide initial MO coefficients for the COLUMBUS calculation. This usually provides a drastic speedup for the CASSCF calculations and helps to ensure that all calculations are based on the same active space. In simple cases it might be acceptable to not provide an initial wavefunction.
Memory
For COLUMBUS, the available memory must be given here.
7.2.5 Input for MOLCAS
Path to MOLCAS
The script will first look for a shell variable $MOLCAS. The user can either confirm that $MOLCAS is the correct executable, or enter the correct path. Note that the setup script will not expand the user ( ) and shell variables (since possibly the ab initio calculations are running on a different machine than the one used for setup). The SHARC– MOLCAS interface will expand and shell variables.
Scratch directory
The script takes the string without any checking. Each individual initial condition uses a subdirectory ICOND_%05i with the appropriate number. and shell variables will only be expanded by the interface during the actual calculation.
Note that you can use, e.g., ./temp/ as scratch directory, which is a subdirectory of the working directory of the interface. Using ./ as scratch directory is not recommended, since the interface will delete the scratch directory.
Template file
Enter the filename for the MOLCAS template input. This file contains the details of the ab initio calculation, like basis set, active orbitals and electrons, number of electrons and number of states for stateaveraging. For details, see the section about the SHARC– MOLCAS interface (6.3). The setup script will check whether the template file contains the necessary entries.
Initial wavefunction file
You can provide initial MO coefficients to MOLCAS. This usually provides a drastic speedup for the CASSCF calculations and helps in ensuring that all calculations are based on the same active space. In simple cases it might be acceptable to not provide an initial wavefunction. For MOLCAS, you have to specify one file for each multiplicity.
7.2.6 Input for analytical potentials
Template file
Enter the filename for the template input specifying the analytical potentials. For details, see the section about the SHARCAnalytical interface (6.5). The setup script will check whether the template file is valid.
7.2.7 Input for Run Scripts
Run script mode
The script setup_init.py generates a run script (Bash) for each initial condition calculation. Due to the large variety of cluster architectures, these run scripts might not work in every case. It is the user’s responsibility to adapt the generated run scripts to his needs.
setup_init.py can generate run scripts for two different schemes how to execute the calculations. With the first scheme, the ab initio calculations are performed in the directory where they were setup (subdirectories of the directory where setup_init.py was started). Note that the interfaces will still use their scratch directories to perform the actual quantum chemistry calculations.
With the second option, the run scripts will transfer the input files for each ab initio calculation to a temporary directory, where the interface is started. After the interface finishes all calculations, the results files are transferred back to the primary directory and the temporary directory is deleted. Note that setup_init.py in any case creates the directory structure in the directory where it was started. The name of the temporary directory can contain shell variables, which will be expanded when the script is running (on the compute host).
Submission script
The setup script can also create a Bash script for the submission of all ab initio calculations to a queueing system. The user has to provide a submission command for that, including any options which might be necessary. This submission script might not work with all queueing systems.
Project name
The user can enter a project name. This is used currently only for the job names of submitted jobs (N option for queueing system).
7.2.8 Output
setup_init.py will create for each initial condition in the given range a directory whose names follow the format ICOND_%05i, where %05i is the index of the initial condition padded with zeroes to 5 digits. Additionally, the directory ICOND_00000 is created for the calculation of the excitation energies at the equilibrium geometry.
To each directory, the following files will be added:
 QM.in: Main input file for the interface, contains the geometry and the control keywords (to specify which quantities need to be calculated).
 run.sh: Run script, which can be started interactively in order to perform the ab initio calculation in this directory. Can also be adapted to a batch script for submission to a queue
 Interfacespecific files: Usually a template file, a file containing additional input for the interface, and an initial wavefunction.
The calculations in each directory can be simply executed by starting run.sh in each directory. In order to perform this task consecutively on a single machine, the script all_run.sh can be executed. The file DONE contains the progress of this calculation.
Alternatively, each run script can be sent to a queueing system (you might need to adapt this script to you cluster system).
In figure 7.1, the directory tree structure setup by setup_init.py is given.
After all calculations are finished, excite.py can be used to collect the results.
7.3 Excitation Selection: excite.py
excite.py has two tasks: adding excitedstate information to the initconds file, and deciding which excited state for which initial condition is a valid initial state for the dynamics.
7.3.1 Usage
The script is interactive, and can be started by simply typing
user@host> python $SHARC/excite.py
7.3.2 Input
Initial condition file
Enter the path to the initial conditions file, to which excite.py will add excitedstate information. This file can already contain excitedstate information (in this case this information can be reused).
Generate excited state list
There are three possibilities to add excitedstate information to the initconds file:
 generate a list of dummy excited states,
 read excitedstate information from the output of the initial ab initio calculations (prepare the calculations with setup_init.py),
 keep the existing excitedstate information in the initconds file.
The first option is mainly used if no initial ab initio calculations need to be performed (e.g., the initial state is known).
In order to use the second option, one should first setup initial excitedstate calculations using setup_init.py (see 7.2) and run the calculations. excite.py can then read the output of the initial calculations and calculate excitation energies and oscillator strengths.
The third option can be used to reuse the information in the initconds file, e.g., to apply a different selection scheme to the states or to just read the number of states.
Path to ab initio results
If excite.py will read the excitedstate information from the ab initio calculation results, here the user has to provide the path to the directory containing the ICOND_%05i subdirectories.
Number of states
If a dummy list of states will be generated, the user has to provide the number of states per multiplicity. Note that a singlet ground state has to be counted as well, e.g. if 4 singlet states are specified, the calculation will involve the S_{0}, S_{1}, S_{2} and S_{3}. Also states of higher multiplicity can be given, e.g. doublet or triplet states (e.g., 2 2 1 for two singlets, two doublets and one triplet).
If the ab initio results are read the number of states will be automatically determined from the results.
Excitedstate representation
When generating new lists of excited states (either dummy states or from ab initio results), the user has to specify the representation of the excited states (either MCH or diagonal representation). The MCH representation is spinfree, meaning that transition dipole moments are only allowed between states of the same multiplicity. For molecules without heavy atoms, this option is sufficient. For heavier atoms, the diagonal representation can be used, which includes the effects of spinorbit coupling on the excitation energies and oscillator strengths.
When reading ab initio results, excite.py will diagonalize the Hamiltonian and transform the transition dipole matrices for each initial condition to obtain the diagonal representation.
When a dummy state list is generated, the representation will only be written to the output file initconds.excited (but has no actual numeric effect for excite.py). Note that the representation which is declared in the initconds.excited file influences how SHARC determines the initial coefficients (see the paragraph on initial coefficients in 4.1.3).
Note that the representation cannot be changed if existing excitedstate information is kept.
Hint: If the ICOND_%05i directories need to be deleted (e.g., due to disk space restrictions), making one readout with excite.py for each representation and saving the results to two different files will preserve all necessary information.
Ionization probabilities
If excite.py detects that the ab initio results contain ionization probabilities, then those can be used instead of the transition dipole moments. Note that in this case the transition dipole moments are not written to the initconds.excited file.
Reference energy
excite.py can read the reference energy (ground state equilibrium energy) directly from the ab initio results. If the ab initio data is read anyways, excite.py already knows the relevant path. If a dummy list of states is generated, the user can provide just the path to the QM.out file of the ab initio calculation for the equilibrium geometry. Otherwise, excite.py will prompt the user to enter a reference energy manually (in hartree).
Initial state selection
Every excited state of each initial condition has a flag specifying it either as a valid initial state or not. excite.py has four modes how to flag the excited states:
 Unselect all excited states,
 User provides a list of initial states,
 States are selected stochastically based on excitation energies and oscillator strengths,
 Keep all existing flags.
The first option can be used if excite.py is only used to read the ab initio results for the generation of an absorption spectrum (using spectrum.py).
The second option can be used to directly specify a list of initial states, if the initial state is known (e.g., starting in the ground state and exciting with an explicit laser field). In this case, the given states of all initial conditions are flagged as initial states.
The third option is only available if excitedstate information exists (i.e., if no dummy list is generated). For details on the stochastic selection procedure, see section 8.5.
The fourth option can only be used if the existing state information is kept. In this case excite.py does nothing except counting the number of flagged initial states.
Excitation window
This option allows to exclude excited states from the selection procedure if they are outside a given energy window. This option is only available if excited state information exists, but not if a dummy list of states is generated (because the dummy states have no defined excitation energy).
For the stochastic selection procedure, states outside the excitation window do not count for the determination of p_{max} (see equation ). This allows to excite, e.g., to a dark nπ^{*} state despite the presence of a much brighter ππ^{*} state.
For the keepflags option, this option can be used to count the number of excited states in the energy window.
Considered states
Here the user can specify the list of desired initial states.
For the stochastic selection procedure, the user can instead exclude certain states from the procedure. Excluded states do not count for the determination of p_{max} (see equation ).
If the number of states per multiplicity is known, excite.py will print a table giving for each state index the multiplicity, quantum number and M_{s} value.
Random number generator seed
The random number generator in excite.py is used in the stochastic selection procedure. Instead of typing an integer, typing “!” will initialize the RNG from the system time. Note that this will not be reproducible, i.e. repeating the excite.py run with “!” as random seed will give a different selection in each run.
7.3.3 Matrix diagonalization
When using the diagonal representation, excite.py needs to diagonalize and multiply matrices. By default, the Python package NUMPY is used, if available. If the script does not find a NUMPY installation, it will use a small Fortran code which comes with the SHARC suite. In order for this to work, you need to set the environment variable $SHARC to the bin/ directory within your SHARC installation. See section 7.13 for more details.
7.3.4 Output
excite.py writes all output to a file <BASE>.excited, where <BASE> is the name of the initial conditions file used as input. The output file is also an initial conditions file, but contains additional information regarding the excited states, the reference energy and the representation of the excited states. An initial conditions file with excitedstate information is needed for the final preparatory step: setting up the dynamics with setup_traj.py.
Additionally, spectrum.py can calculate absorption spectra from excitedstate initial condition files.
7.3.5 Specification of the initconds.excited file format
The initial conditions files initconds and initconds.excited contain lists of initial conditions, which are needed for the setup of trajectories. An initial condition is a set of initial coordinates of all atoms and corresponding initial velocities of each atom, and optionally a list of excited state informations. In the following, the format of this file is specified.
The file contains of a header, followed by the body of the file containing a list of the initial conditions.
File header
An examplary header looks like:
SHARC Initial conditions file, version 0.2 <Excited> Ninit 100 Natom 2 Repr MCH Eref 0.50 Eharm 0.04 States 2 0 1 Equilibrium H 1.0 0.0 0.0 0.0 1.00782503 0.0 0.0 0.0 H 1.0 1.5 0.0 0.0 1.00782503 0.0 0.0 0.0
The first line must read SHARC Initial conditions file, version <VERSION>, with the correct version string followed. The string Excited is optional, and marks an initial conditions file as being an output file of excite.py (setup_traj.py will only accept files marked like this). The following lines contain:
 the number of initial conditions,
 the number of atoms,
 the electronic state representation (a string which is None, MCH or diag),
 the reference energy (hartree),
 the harmonic energy (zero point energy in the harmonic approximation, hartree),
 optionally the number of states per multiplicity.
After the header, first the equilibrium geometry is expected. It is demarked with the keyword Equilibrium, followed by n_{atom} lines, each specifying one atom. Unlike the actual initial conditions, the equilibrium geometry does not have a list of excited states or defined energies.
File body
The file body contains a list of initial conditions. Each initial condition is specified by a block starting with a line containing the string Index and the number of the initial condition. In the file, the initial conditions are expected to appear in order.
A block specifying an initial condition looks like:
Index 1 Atoms H 1.0 0.02 0.0 0.0 1.00782503 0.001 0.0 0.0 H 1.0 1.52 0.0 0.0 1.00782503 0.001 0.0 0.0 States 001 0.49 0.49 0.16 0.0 0.03 0.0 0.05 0.0 0.0 0.00 False 002 0.25 0.49 0.02 0.0 0.43 0.0 1.77 0.0 6.5 0.53 True 003 0.40 0.49 0.00 0.0 0.00 0.0 0.00 0.0 2.5 0.00 False 004 0.40 0.49 0.00 0.0 0.00 0.0 0.00 0.0 2.5 0.00 False 005 0.40 0.49 0.00 0.0 0.00 0.0 0.00 0.0 2.5 0.00 False Ekin 0.004 a.u. Epot_harm 0.026 a.u. Epot 0.013 a.u. Etot_harm 0.030 a.u. Etot 0.018 a.u.
The formal structure of such a block is as follows. After the line containing the keyword Index and the index number, the keyword Atoms indicates the start of the list of atoms. Each atom is specified on one line:
 symbol,
 nuclear charge,
 x, y, z coordinate in Bohrs,
 atomic mass,
 x, y and z component of nuclear velocity in atomic units.
After the atom list, the keyword States indicates the list of electronic states. This list consists of one line per electronic state, but can be empty, if no information of the electronic states is available. Each line consists of:
 state number (starting with 1),
 state energy in Hartree,
 reference energy in Hartree (usually the energy of the lowest state),
 six numbers defining the transition dipole moment to the reference state (usually the lowest state),
 the excitation energy in eV,
 the oscillator strength,
 a string which is either True or False, specifying whether the electronic state was selected by excite.py as initial electronic state.
The transition dipole moments are specified by six floating point numbers, which are real part of the x component, imaginary part of the x component, then the real and imaginary parts for the y and finally the z component (the transition dipole moments can be complex in the diagonal representation).
The electronic state list is terminated with the keyword Ekin, which at the same time gives the kinetic energy of all atoms. The remaining entries give the potential energy in the harmonic approximation and the actual potential energy, as well as the total energy.
7.4 Setup of Trajectories: setup_traj.py
This interactive script prepares the input for the excitedstate dynamics simulations with SHARC. It works similarly to setup_init.py, reading an initial conditions file, prompting the user for a number of input parameters, and finally prepares one directory per trajectory. However, the setup_traj.py input section is noticeably longer, because many options for the dynamics are covered.
7.4.1 Input
Initial conditions file
Please be aware that setup_traj.py needs an initial conditions file generated by excite.py (Files generated by wigner.py are not allowed). The script reads the number of initial states, the representation and the reference energy automatically from the file.
Number of states
This is the total number of states per multiplicity included in the dynamics calculation. Affects the keyword nstates in the SHARC input file.
Only advanced users should use here a different number of states than given to setup_init.py. In this case, the excitedstate information in the initial conditions file might be inconsistent. For example, if 10 singlets and 10 triplets were included in the initial calculations, but only 5 singlets and 5 triplets in the dynamics, then the sixth entry in the initial conditions file corresponds to S_{5}, while setup_traj.py assumes the sixth entry to correspond to T_{1}.
Active states
States can be frozen for the dynamics calculation here. See section 8.2 for a general description of state freezing in SHARC. Only the highest states in each multiplicity can be frozen, it is not possible to, e.g., freeze the ground state in simulations where ground state relaxation is negligible. Affects the keyword actstates.
Contents of the initial conditions file
Optionally, a map of the contents of the initial conditions file can be displayed during the execution of setup_traj.py, showing for each state which initial conditions were selected (and which initial conditions do not have the necessary excitedstate information). For each state, a table is given, where each symbol represents one initial condition. A dot “.” represents an initial condition where information about the current excited state is available, but which is not selected for dynamics. A hash mark “#” represents an initial condition which is selected for dynamics. A question mark “?” represents initial conditions for which no information about the excited state is available (e.g. if the initial excitedstate calculation failed). The tutorial shows an example of this output.
The content of the initial conditions file is also summarized in a table giving the number of initial conditions selected per state.
Initial states for dynamics setup
The user has to input all states from which trajectories should be launched. The numbers must be entered according to the above table giving the number of selected initial conditions per state. It is not allowed to specify inactive states as initial states. The script will give the number of trajectories which can be setup with the specified set of states. If no trajectories can be setup, the user has to specify another set of initial states. The initial state will be written to the SHARC input, specified in the same representation as given in the initial conditions file. The initial coefficients will be determined automatically by SHARC, according to the description in section 4.1.3.
Starting index for dynamics setup
Specifies the first initial condition within the initial condition file to be included in the setup. This is useful, for example, if the user might setup 50 trajectories starting with index 1. setup_traj.py reports afterwards the last initial condition to be used for setup, e.g. index 90. Later, the user can setup additional trajectories, starting with index 91.
Random number generator seed
The random number generator in setup_traj.py is used to randomly generate RNG seeds for the SHARC input. Instead of typing an integer, typing “!” will initialize the RNG from the system time. Note that this will not be reproducible, i.e. repeating the setup_traj.py run (with the same input) with “!” as random seed will give for the same trajectories different RNG seeds. Affects the keyword RNGseed.
Interface
In this point, choose any of the displayed interfaces to carry out the ab initio calculations. Enter the corresponding number. The choice of the interface influences some dynamics options which can be set in the next section of the setup_traj.py input.
Simulation Time
This is the maximum time that SHARC will run the dynamics simulation. If trajectories need to be run for longer time, it is recommended to first let the simulation finish. Afterwards, increase the simulation time in the corresponding SHARC input file (keyword tmax) and add the restart keyword (also make sure that the norestart keyword is not present). Then the simulation can be restarted by running again the run.sh script. Sets the keyword tmax in the SHARC input files.
Simulation Timestep
This gives the timestep for the dynamics. The onthefly ab initio calculations are performed with this timestep, as is the propagation of the nuclear coordinates. A shorter timestep gives more accurate results, especially if light atoms (hydrogen) are subjected to high kinetic energies or steep gradients. Of course a shorter timestep is computationally more expensive. A good compromise in many situations is 0.5 fs. Sets the keyword stepsize in the SHARC input files.
Number of substeps
This gives the number of substeps for the interpolation of the Hamiltonian for the propagation of the electronic wavefunction. Usually, 25 substeps are sufficient. In cases where the diagonal elements of the Hamiltonian are very large (very large excitation energies or a badly chosen reference energy) more substeps are necessary. Sets the keyword nsubsteps in the SHARC input files.
Prematurely terminate trajectories
Usually, trajectories which relaxed to the ground state do not recross to an excited state, but vibrate indefinitely in the ground state. If the user is not interested in these vibrations, such trajectories can be terminated prematurely in order to save computational ressources. A threshold of 1020 fs is usually a good choice to safely detect ground state relaxation. Sets the keyword killafter in the SHARC input files.
Representation for the dynamics
Either the diagonal representation can be chosen (by typing “yes“) to perform dynamics with the SHARC methodology, or the dynamics can be performed on the MCH states (spindiabatic dynamics [24], FISH [23]). Sets the keyword surf in the SHARC input files.
Nonadiabatic couplings
Electronic propagation can be performed with temporal derivatives, nonadiabatic coupling vectors or overlap matrices (Local diabatization). Enter the corresponding number. Note that depending on the chosen interface, some options might not be available, as displayed by setup_traj.py. Sets the keyword coupling in the SHARC input files.
Gradient transformation
The nonadiabatic coupling vectors can be used to correctly transform the gradients to the diagonal representation. If nonadiabatic coupling vectors are used anyways, this option is strongly recommended, since it gives more accurate gradients for no additional cost. Sets the keyword gradcorrect in the SHARC input files. If the dynamics uses the MCH representation, this question is not asked.
Surface hop treatment
This option determines how the total energy is conserved after a surface hop. Sets the keyword ekincorrect in the SHARC input files.
Decoherence correction
For most applications, a decoherence correction should be enabled. Requesting decoherence correction here leads to the addition of the keyword decoherence to the SHARC input files.
Note that setup_traj.py does not allow to modify the α value. In order to change the keyword decoherence_param, the user has to manually edit the SHARC input files.
Scaling and Damping
These two prompts set the keywords scaling and damping in the SHARC input files. The scaling parameter has to be positive, and the damping parameter has to be in the interval [0,1].
Gradient and nonadiabatic coupling selection
For dynamics in the MCH representation, selection of gradients is used by default, and only one gradient (of the current state) is calculated. Selection of nonadiabatic couplings is only relevant if they are used (for propagation, gradient correction or rescaling of the velocities after a surface hop). For the selection threshold, usually 0.5 eV is sufficient, except if spinorbit coupling is very strong and hence the gradients mix strongly.
Sets the keywords grad_select and nac_select in the SHARC input files.
Laser file
The user can specify to use an external laser field during the dynamics, and has to provide the path to the laser file (see section 7.5 and 4.5). setup_traj.py will check whether the number of steps and the timesteps are compatible to the dynamics. If the interface can provide dipole moment gradients, setup_traj.py will also ask whether dipole moment gradients should be included in the simulations.
7.4.2 Interfacespecific input
This input section is basically the same as for setup_init.py (sections 7.2.3 to 7.2.6). Note that for the dynamics simulations an initial wavefunction file is even more strongly recommended than for the initial excitedstate calculations.
For MRCI dynamics using the COLUMBUS interface, for performance reasons it is strongly advisable to obtain the excitlistfiles before starting the dynamics.
7.4.3 Run script setup
Also this input section is very similar to the one in setup_init.py (see section 7.2).
7.4.4 Output
setup_traj.py will create for each initial state a directory where all trajectories starting in this state will be put. If the initial conditions file specified that the initial conditions are in the MCH representation, then the initial states will be assumed to be in the MCH representation as well. In this case, the directories will be named Singlet_0, Singlet_1, …, Doublet_0, Triplet_1, … If the initial states are in the diagonal representation, then the directories are simply called X_1, … since they do not have a definite spin.
In each directory, subdirectories called TRAJ_%05i are created, where %05i is the initial condition index, padded to 5 digits with zeroes. In each trajectory’s directory, an SHARC input file called input will be created, which contains all the dynamics options chosen during the setup_traj.py run. Also, files geom and veloc will be created. For trajectories setup with setup_traj.py, the determination of the initial wavefunction coefficients is done by SHARC.
Furthermore, in each trajectory directory a subdirectory QM is created, where the runQM.sh script containing the call to the interface is put. In the directory QM also all interfacespecific input files will be copied.
For each trajectory, a run.sh script will be created, which can be executed to run the dynamics simulation. You might need to adapt the run script to your cluster setup.
Optionally, setup_traj.py also creates a script all_qsub.sh, which can be executed to submit all trajectories to a queueing system. You might need to adapt also this script to your cluster setup.
The full directory structure created by setup_traj.py is given in figure 7.2.
7.5 Laser field generation: laser.x
The Fortran code laser.x can generate files containing laser fields which can be used with SHARC. It is possible to superimpose several lasers, use different polarizations and apply a number of chirp parameters.
7.5.1 Usage
The program is simply called by
user@host> python $SHARC/laser.x
It will interactively ask for the laser parameters. After input is complete, it writes the laser field to the file laser in the format which SHARC expects (see 4.5).
Similar to the interactive Python scripts, laser.x will also write the user input to KEYSTROKES.laser. After modifying this file, it can be used to directly execute laser.x without doing the interactive input again:
user@host> python $SHARC/laser.x < KEYSTROKES.laser
7.5.2 Input
The first four options are global and need to be entered only once, all remaining input options need to be given for every laser pulse. For the definition of laser fields see section 8.8.
Number of lasers
Any number of lasers can be used. The output file will contain the sum of all laser pulses defined.
Realvalued field
If this is true, the output file will only contain the real parts of the laser field, while the columns defining the imaginary part of the field will be zero. Note, however, that SHARC will anyways only use the real part of the field in the simulations.
Time interval and steps
The definitions of the starting time, end time and time step of the laser field must exactly match the simulation time and time substeps of the SHARC simulation. Note, that the laser field must always start at t=0 fs to be used with SHARC. The end time for the laser field must therefore coincide with the total simulation time given in the SHARC input. The number of time steps for the laser field is t_{total}/∆t_{sub} +1.
Files for debugging
This option is normally not needed, and can be set to False. If set to True, the chirped and unchirped laser fields in both time and frequency domain will be written to files called DEBUG_... .
Polarization vector
The polarization vector p (will be normalized).
Type of envelope
There are two options possible for the envelope function E(t), either a Gaussian envelope or a sinusoidal one (see 8.8).
Field strength
There are two input lines for the field strength E_{0}, the first defining the unit in which the field strength is defined, the second gives the corresponding number. Field strength can be read in in GV/m, TW/cm^{−2} or atomic units.
FWHM and time intervals
This option depends on the type of envelope chosen. While in both cases all 5 numbers need to be entered, for a Gaussian pulse only the first and third number have an effect. For a sinusoidal pulse all but the first number has an effect.
For a Gaussian pulse, the first argument corresponds to FWHM in equation and the third argument to t_{c} in .
For a sinusoidal pulse, the second, third, fourth and fifth argument correspond to t_{0}, t_{c}, t_{c2} and t_{e}, respectively, in equation .
Central frequency
There are two input lines for the central frequency ω_{0}. The first defines the unit (wavelength in nm, energy in eV, or atomic units). The second line gives the value.
Phase
The total phase ϕ is given in multiples of π. For example, the input “1.5” gives a phase of [(3π)/2].
Chirp parameters
There are four lines giving the chirp parameters b_{1}, b_{2}, b_{3} and b_{4}. See equation for the meaning of these parameters.
7.6 Calculation of Absorption Spectra: spectrum.py
Aside from setting up trajectories, the initconds.excited files can also be used to generate absorption spectra based on the excitation energies and oscillator strengths in the file. The script spectrum.py calculates Gaussian or Lorentzian convolutions of these data in order to obtain spectra. See section 8.1 for further details.
spectrum.py evaluates the absorption spectrum on a grid for all states it finds in an initial conditions file. Using commandline options, some initial conditions can be omitted in the convolution, see table 7.2.
7.6.1 Input
The script is executed with the initial conditions file as argument:
user@host> python $SHARC/spectrum.py [OPTIONS] initconds.excited
The script accepts a number of commandline options, which are given in table 7.2.
Option  Description  Default 
h  Display help message and quit.  – 
o FILENAME  Output filename (for the spectrum)  spectrum.out 
n INTEGER  Number of grid points  500 
e FLOAT FLOAT  Energy range (eV) for the spectrum  1 to 10 eV 
i INTEGER INTEGER  Index range for the initial conditions  1 to 1000 
f FLOAT  FWHM (eV) for the spectrum  0.1 eV 
G  Gaussian convolution  Gaussian 
L  Lorentzian convolution  Gaussian 
s  Use only selected initial conditions  Use all 
l  Make a line spectrum  Convolution 
–gnuplot FILENAME  Write a GNUPLOT script  No GNUPLOT script 
7.6.2 Output
The script writes the absorption spectrum to a file (by default spectrum.out). Using the o option, the user can redirect the output to a suitable file. The output is a table containing n+2 columns, where n is the number of states found in the initial conditions file. The first column gives the energy in eV, within the given energy interval. In columns 2 to n+1 the statewise absorption spectra are given. The last column contains the total absorption spectrum, i.e., the sum over all states. The table has n_{grid}+1 rows. For line spectra the output format is exactly the same, however, the file will contain one row for each excited state of each initial condition in the initial conditions file.
Additionally, the script writes some information about the calculation to standard output, among these the maximum of the spectrum, which can be used in order to normalize the spectrum. The reported maximum is simply the largest value in the last column of the spectrum.
If requested, the script generates a GNUPLOT script, which can be used to directly plot the spectrum.
7.7 File transfer: retrieve.sh
Usually, SHARC will run on some temporary directory, and not in the directory where the trajectories have been submitted from. The shell script retrieve.sh is a simple scp wrapper, which can be executed (in a directory where a trajectory has been sent from) in order to retrieve the output files of this trajectory. This might not work for every cluster setup.
It relies on the presence of the file host_infos. All trajectories set up with setup_traj.py create this file after the trajectory has been started with run.sh. retrieve.sh reads host_infos to determine the hostname and working directory of the trajectory and then uses scp to retrieve the output and restart files.
The script can be called with the option “lis” in order to only retrieve the output.lis file, but not the other output files.
If the script is called with the option “res” then also the restart files and the content of the restart/ directory are copied.
It is advisable to configure publickey authentification for the hosts running the trajectories, so that not for every execution of retrieve.sh a password has to be entered.
7.8 Data Extractor: data_extractor.x
The output of SHARC is mainly written to output.dat. In order to obtain plottable files in tabular format, the Fortran program data_extractor.x is used.
7.8.1 Usage
The data_extractor.x is a command line tool, and is called with the output.dat file as an argument.
user@host> $SHARC/data_extractor.x output.dat
The program will create a directory output_data/ in the current working directory (not necessarily in the directory where output.dat resides). In this directory, several files are written, containing e.g. the potential energies depending on time, etc.
The program will extract the complete output.dat file until it reaches the EOF. No further input except for the data file can be given.
7.8.2 Output
After the program finishes, the directory output_data/ contains a number of files. In each file, the number of columns is dependent in the total number n of states i ∈ {1…n}. The content of the files is listed in Table 7.3.
The file expec.out contains the information of energy.out, spin.out and fosc.out in one file. The content of expec.out can be conveniently plotted by using make_gnuscript.py to generate a GNUPLOT script.
File  # Columns  Columns  
coeff_diag.out  2+2n  1  Time t (fs) 
2  Norm of wavefunction ∑_{j} c_{j}^{diag}^{2}  
1+2j  ℜ(c_{j}^{diag})  
2+2j  ℑ(c_{j}^{diag})  
coeff_MCH.out  2+2n  1  Time t (fs) 
2  Norm of wavefunction ∑_{j} c_{j}^{MCH}^{2}  
1+2j  ℜ(c_{j}^{MCH})  
2+2j  ℑ(c_{j}^{MCH})  
energy.out  4+n  1  Time t (fs) 
2  Kinetic energy (eV)  
3  Potential energy (eV) of active state  
4  Total energy (eV)  
4+j  Potential energy (eV) of state j  
fosc.out  2+n  1  Time t (fs) 
2  Oscillator strength of active state  
2+j  Oscillator strength of state j  
spin.out  2+n  1  Time t (fs) 
2  Total spin expectation value of active state  
2+j  Total spin expectation value of state j  
prob.out  2+n  1  Time t (fs) 
2  Random number from surface hopping  
2+j  Cumulated hopping probability ∑_{k=1}^{j} P_{k}  
expec.out  4+3n  1  Time t (fs) 
2  Kinetic energy (eV)  
3  Potential energy (eV) of active state  
4  Total energy (eV)  
4+j  Potential energy (eV) of state j  
4+n+j  Total spin expectation value of state j  
4+2n+j  Oscillator strength of state j  
7.9 Plotting the Extracted Data: make_gnuscript.py
The contents of the output files of data_extractor.x can be plotted with GNUPLOT. In order to quickly generate an appropriate GNUPLOT script, make_gnuscript.py can be used. The usage is:
user@host> python $SHARC/make_gnuscript.py <S> [<D> [<T> [<Q> ... ] ] ]
make_gnuscript.py takes between 1 and 8 integers as commandline arguments, specifying the number of singlet, doublet, triplet, etc. states. It writes an appropriate GNUPLOT script to standard out, hence redirect the output to a file, e.g.:
user@host> python $SHARC/make_gnuscript.py 3 0 2 > gnuscript.gp
Then, GNUPLOT can be run in the output_data directory of a trajectory:
user@host> gnuplot gnuscript.gp
This can also be accomplished in one step using a pipe, e.g.:
user@host> python $SHARC/make_gnuscript.py 3 0 2  gnuplot
The created plot script generates four different plots (press ENTER in the commandline where you started GNUPLOT to go to the next plot). The first plot shows the potential energy of all states in the dynamics over time in the diagonal representation. The currently occupied state is marked with black circles. A thin black line gives the total energy (sum of the kinetic energy and the potential energy of the currently occupied state). Each state is colored, with one color as contour and one color at the core of the line. The contour color represents the total spin expectation value of the state. The core color represents the oscillator strength of the state with the lowest state. See figure 7.3 for the relevant color code.
Note that by definition the “oscillator strength” of the lowest state with itself is exactly zero, hence the lowest state is also light grey. This dual coloring allows for a quick recognition of different types of states in the dynamics, e.g. singlets vs. triplets or nπ^{*} vs. ππ^{*} states.
The second plot shows the population c_{i}^{MCH}^{2} of the MCH electronic states over time. The line colors are autogenerated in order to give a large spread of all colors over the excited states, but the colors might be suboptimal, e.g. for printing. In this cases, the user should manually adjust the colors in the generated script.
The third plot shows the population c_{i}^{diag}^{2} of the diagonal electronic states over time. These are the populations which are actually used for surface hopping. However, since these states are spinmixed, it is usually difficult to interpret these populations.
The fourth plot shows the surface hopping probabilities over time. The plot is setup in such a way that the visible area corresponding to a certain state is proportional to the probability to hop into the state. Hence, if for a given timestep the random number (black circles) lies within a colored area, a surface hop to the corresponding state is performed.
7.10 Calculation of Ensemble Populations: populations.py
For an ensemble of trajectories, usually one of the most relevant results are ensembleaveraged populations. The interactive script populations.py collects these populations from a set of trajectories.
Different methods to obtain populations or quantities approximating populations can be collected, as described below.
7.10.1 Usage
The script is interactive, simply start it with no commandline arguments or options:
user@host> python $SHARC/populations.py
Depending on the analysis mode (see below) it might be necessary to run data_extractor.x for each trajectory prior to running populations.py.
Paths to trajectories
First the script asks the user to specify all directories for whose content the analysis should be performed. Enter one directory path at a time, and finish the directory input section by typing “end“. Please do not specify each trajectory directory separately, but specify their parent directories, e.g. the directories Singlet_1 and Singlet_2. populations.py will automatically include all trajectories contained in these directories.
If you want to exclude certain trajectories from the analysis, it is sufficient to create an empty file called CRASHED or RUNNING in the corresponding trajectory directory. populations.py will ignore all directories containing one of these files.
Analysis mode
Using populations.py, there are two basic ways in obtaining the excitedstate populations. The first way is to count the number of trajectories for which a certain condition holds. For example, the number of trajectories in each classical state can be obtained in this way. However, it is also possible to count the number of trajectories for which the total spin expectation value is within a certain interval.
The second way to obtain populations is to obtain the sum of the absolute squares of the quantum amplitudes over all trajectories. Table 7.4 contains a list of all possible analysis modes.
Mode  Description  From which  Extract? 
file?  
1  For each diagonal state count how many trajectories have this state as active state.  output.lis  No 
2  For each MCH state count how many trajectories have this state as approximate active state (see section 8.11.1).  output.lis  No 
3  For each MCH state count how many trajectories have this state as approximate active state (see section 8.11.1). Multiplet components are summed up.  output.lis  No 
4  Generate a histogram with definable bins (variable width). Bin the trajectories according to their total spin expectation value (of the currently active diagonal state).  output.lis  No 
5  Generate a histogram with definable bins (variable width). Bin the trajectories according to their state dipole moment expectation value (of the currently active diagonal state).  output.lis  No 
6  Generate a histogram with definable bins (variable width). Bin the trajectories according to the oscillator strength between lowest and currently active diagonal states.  output_data/fosc.out  Yes 
7  Calculate the sum of the absolute squares of the diagonal coefficients for each state.  output_data/coeff_diag.out  Yes 
8  Calculate the sum of the absolute squares of the MCH coefficients for each state.  output_data/coeff_MCH.out  Yes 
9  Calculate the sum of the absolute squares of the MCH coefficients for each state. Multiplet components are summed up.  output_data/coeff_MCH.out  Yes 
10  Calculate the sum of the absolute squares of the diabatic coefficients for each state (Only for trajectories with local diabatization).  output_data/coeff_diab.out  Yes 
Run data extractor
For analysis modes 6, 7, 8, 9 and 10 it is necessary to first run the data extractor (see section 7.8). This task can be accomplished by populations.py. However, for a large ensemble or for long trajectories this may take some time. Hence, it is not necessary to perform this step each time populations.py is run.
populations.py will detect whether the file output.dat or the content of output_data/ is more recent. Only if output.dat is newer the data_extractor.x will be run for this trajectory.
Note that mode 10 can only be used for trajectories using local diabatization propagation (keyword coupling overlap in SHARC input file).
Number of states
For analysis modes 1, 2, 3, 7, 8 and 9 it is necessary to specify the number of states in each multiplicity. The number is autodetected from the input file of one of the trajectories.
Intervals
For analysis modes 4, 5 and 6 the user must specify the intervals for the classification of the trajectories. The user has to input a list of interval borders, e.g.:
Please enter the interval limits, all on one line. Interval limits: 1e3 0.01 0.1 1
Note that scientific notation can be used. Based on this input, for each timestep a histogram is created with the number of trajectories in each interval. The histogram bins are:
 x ≤ 10^{−3}
 10^{−3} < x ≤ 10^{−2}
 10^{−2} < x ≤ 10^{−1}
 10^{−2} < x ≤ 10^{0}
 10^{0} < x
Note that there is always one more bin that interval borders entered.
Normalization
If desired, populations.py can normalize the populations by dividing the populations by the number of trajectories.
Maximum simulation time
This gives the maximum simulation time until which the populations are analyzed. For trajectories which are shorter than this value, the last population information is used to make the trajectory long enough. Trajectories which are longer are not analyzed to the end. populations.py prints the length of the shortest and longest trajectories after the analysis.
Gnuplot script
populations.py can generate an appropriate GNUPLOT script for the performed population analysis.
7.10.2 Output
By default, populations.py writes the resulting populations to pop.out. If the file already exists, the user is ask whether it shall be overwritten, or to provide an alternative filename. Note that the output file is checked only after the analysis is completed, so the program might run for a considerable amount of time before asking for the output file.
7.11 Obtaining special geometries: crossing.py
In many cases, it is also important to obtain certain special geometries from the trajectories. The script crossing.py extracts geometries fulfilling special conditions from an ensemble of trajectories.
Currently, crossing.py finds geometries where the approximate MCH state (see section 8.11.1) of the last timestep is different from the MCH state of the current timestep (i.e. crossing.py finds geometries where surface hops occured).
7.11.1 Usage
The script is interactive, simply start it with no commandline arguments or options:
user@host> python $SHARC/crossing.py
The input to the script is very similar to the one of populations.py.
Paths to trajectories
First the script asks the user to specify all directories for whose content the analysis should be performed. Enter one directory path at a time, and finish the directory input section by typing “end“. Please do not specify each trajectory directory separately, but specify their parent directories, e.g. the directories Singlet_1 and Singlet_2. crossing.py will automatically include all trajectories contained in these directories.
If you want to exclude certain trajectories from the analysis, it is sufficient to create an empty file called CRASHED or RUNNING in the corresponding trajectory directory. crossing.py will ignore all directories containing one of these files.
Analysis mode
Currently, crossing.py only supports one analysis mode, where crossing.py is scanning for each trajectory the file output.lis. If the occupied MCH state (column 4 in output file output.lis) changes from one timestep to the next, it is checked whether the old and new MCH states are the ones specified by the user. If this is the case, the geometry corresponding to the new timestep (t) is retrieved from output.xyz (lines t(n_{atom}+2)+1 to t(n_{atom}+2)+n_{atom}).
States involved in surface hop
First, the user has to specify the permissible old MCH state. The state has to be specified with two integers, the first giving the multiplicity (1=singlet, …), the second the state within the multiplicity (1 1=S_{0}, 1 2=S_{1}, etc.). If a state of higher multiplicity is given, crossing.py will report all geometries where the old MCH state is any of the multiplet components.
For the new MCH state, the same is valid.
Third, the direction of the surface hop has to be specified. Choosing “Backwards” has the same effect as exchanging the old and new MCH states.
7.11.2 Output
All geometries are in the end written to an output file, by default crossing.xyz. The file is in standard xyz format. The comment of each geometry gives the path to the trajectory where this geometry was extracted, the simulation time and the diagonal and MCH states at this simulation time.
7.12 Internal Coordinates Analysis: geo.py
SHARC writes at every timestep the molecular geometry to the file output.xyz. The noninteractive script geo.py can be used in order to extract internal coordinates from xyz files. The usage is:
user@host> python $SHARC/geo.py [options] < Geo.inp > Geo.out
By default, the coordinates are read from output.xyz, but this can be changed with the g option (see table 7.6). Note that the internal coordinate specifications are read from standard input and the result table is written to standard out.
7.12.1 Input
The specifications for the desired internal coordinates are read from standard input. It follows a simple syntax, where each internal coordinate is specified by a single line of input. Each line starts with a oneletter key which specifies the type of internal coordinate (e.g. bond length, angle, dihedral, …). The key is followed by a list of integers, specifying which atoms should be measured. As a simple example, r 1 2 specifies the bond length (r is the key for bond lengths) between atoms 1 and 2. Note that the numbering of the atoms starts with 1. Each line of input is checked for consistency (whether any atom index is larger than the number of atoms, repeated atom indices, misspelled keys, wrong number of atom indices, …), and erroneous lines are ignored (this is indicate by an error message).
Table 7.5 lists the available types of internal coordinates. The output is a table, where the first column is the time (Actually, the geometries are just enumerated starting with zero, and the number multiplied by the timestep from the t option). The successive columns in the output table list the results of the internal coordinates calculations. Each request generates at least one column, see table 7.5.
Note that for most internal coordinates, the order of the atoms is crucial, since e.g. a_{123} ≠ a_{213}. This also holds for the CremerPople parameter requests. For these input lines, the atoms should be listed in the order they appear in the ring (clockwise or counterclockwise).
Key  Atom  Description  Output columns 
Indices  
x  a  x coordinate of atom a  x 
y  a  y coordinate of atom a  y 
z  a  z coordinate of atom a  z 
r  a b  Bond length between a and b  r 
a  a b c  Angle between a–b and b–c  a 
d  a b c d  Dihedral (Angle between planes (a,b,c)  d 
and (b,c,d))  
p  a b c d  Pyramidalization angle (90^{°} minus angle  p 
between bond a–b and plane (b,c,d))  
5  a b c d e  CremerPople parameters for 5membered  q_{2}, ϕ_{2} 
rings [40].  
6  a b c d e f  CremerPople parameters for 6membered  Q, ϕ, θ, Boeyens 
rings [40] and Boeyens classification [41].  
c  Writes the comment (second line of the  Comment  
xyz format) to the table.  
As an advice, it is always a good idea to put the comment as the last request, if at all. Since the comment may contain blanks, having the comment not as the very last column might make it impossible to plot the resulting table.
The Boeyens classification symbols which are output for 6membered rings are reported in L^{A}T_{E}X math code. Note that in the Boeyens classification scheme a number of symbols are equivalent, and only one symbol is reported. For example, by definition the chair configuration ^{1}C_{4} is equivalent to ^{5}C_{2} and ^{3}C_{6}.
7.12.2 Options
geo.py accepts a number of commandline options, see table 7.6. All options have sensible defaults. However, especially if long comments should be written to the output file, it might be necessary to increase the field width. Note that the minimum column width is 20 so that the table header can be printed correctly.
Option  Description  Default 
h  Display help message and quit.  – 
b  Report x, y, z, r, q_{2} and Q in Bohrs  Angstrom 
r  Report a, d, p, ϕ_{2}, ϕ and θ in Radians  Degrees 
g FILENAME  Read coordinates from the specified file  output.xyz 
t FLOAT  Assumed timestep between successive geometries (fs)  1.0 fs 
p INTEGER  Precision (Number of decimals, max. width3)  4 
w INTEGER  Width of each column (min. 20)  20 
7.13 Diagonalization Helper: diagonalizer.x
The small program diagonalizer.x is used by the Python scripts to diagonalize matrices if the NUMPY library is not available. Currently, only excite.py and SHARC_Analytical.py need to diagonalize matrices. The program diagonalizer.x is implemented as a simple frontend to the LAPACK libary.
The program reads from stdin. The first line consists of the letter “r” or “c”, followed by two integers giving the matrix dimensions. For “r”, the program assumes a real symmetric matrix, for “c” a Hermitian matrix. The matrix must be square.
The second line is a comment and is ignored.
In the following lines, the matrix elements are given. On each line one row of the matrix has to be written. If the matrix is complex, each line contains the double number of entries, where real and imaginary part are always given subsequently.
The following is an example input:
c 2 2 comment 1.0 0.0 2.0 0.1 2.0 0.1 6.0 0.0
The diagonal matrix and the matrix of eigenvectors is written to stdout.
Chapter 8 Methodology
In this chapter different aspects of SHARC simulations are discussed in detail. The topics are ordered alphabetically.
8.1 Absorption Spectrum
Using spectrum.py, an absorption spectrum can be calculated as the sum over the absorption spectra of each individual initial condition:
where i runs over the initial conditions.
The spectrum of a single initial condition is the convolution of its line spectrum, defined through a set of tuples (E^{α},f_{osc}^{α}) for each electronic state α, where E^{α} is the excitation energy and f_{osc}^{α}) is the oscillator strength.
The convolution of the line spectrum can be performed with spectrum.py using either Gaussian or Lorentzian functions. The contribution of a state α to the absorption spectrum σ^{α}(E) is given by:
or
where FWHM is the full width at half maximum.
8.2 Active and inactive states
SHARC allows to “freeze” certain states, which then do not participate in the dynamics. Only energies and dipole moments are calculated, but all couplings are disabled. In this way, these states are never visited (hence also no gradients and nonadiabatic couplings are calculated, making the inclusion of these states cheap). Example:
nstates 2 0 2
actstates 2 0 1
In the example given, state T_{2} is frozen. The corresponding Hamiltonian looks like:
where all matrix elements marked
red 
are deleted, since T_{2} is frozen.
The corresponding matrix elements are also deleted from the nonadiabatic coupling and overlap matrices. For propagation including laser fields, also the corresponding transition dipole moments are neglected, while the transition dipole moments still show up in the output (in order to characterize the frozen states).
Active and frozen states are defined with the states and actstates keywords in the input file. Note that only the highest states in each multiplicity can be frozen, i.e., it is not possible to freeze the T_{1} while having T_{2} active. However, it is possible to freeze all states of a certain multiplicity.
8.3 Damping
If damping is activated in SHARC (keyword dampeddyn), in each timestep the following modification to the velocity vector is made
where C is the damping factor given in the input. Hence, in each timestep the kinetic energy is modified by
The damping factor C must be between 0 and 1.
8.4 Decoherence
In surface hopping, without any corrections the coherence between the states is usually too large [19]. A trajectory in state β, but where state α has a large coefficient, will still travel according to the gradient of state β. However, the gradients of state α are almost certainly different to the ones of state β. As a consequence, too much population of state α is following the gradient of state β. Decoherence corrections damp in different ways the population of all states α ≠ β, so that only population of β follows the gradient of state β.
The decoherence correction currently implemented in SHARC is based on the energies of the excited states (“energybased decoherence”, or EDC, as given in [50]). After the surface hopping procedure, when the system is in state β, the coefficients are updated by the following relation
where C is the decoherence parameter. The decoherence correction can be activated with the keyword decoherence in the input file. The decoherence parameter C can be set with the keyword decoherence_param (the default is 0.1 hartree, as suggested in [50]).
8.5 Excitation Selection
excite.py can select initial active states for the dynamics based on for each initial condition k the excitation energies E_{k,α} and the oscillator strengths f^{osc}_{k,α} are known for the excited states α. The procedure is based on reference [39].
First, for all excited states of all initial conditions, the maximum value p_{max} of
is found. Then for each excited state, a random number r_{k,α} is picked from [0,1]. If
then the excited state is selected for performing the dynamics. This excitedstate selection scheme is taken from [39].
Within excite.py it is possible to restrict the selection to a subset of all excited states. In this case, also p_{max} is only determined based on this subset of excited states.
8.6 Internal coordinates definitions
In this section, the internal coordinates available in geo.py are defined.
Bond length
The bond length between two atoms a and b is:
Bond angle
The bond angle for atoms a, b and c is defined as the angle
where v_{ba} is the vector from atom b to atom a.
Dihedral
The dihedral angle is defined via the vectors w_{1} and w_{2}:
The dihedral is given as the angle between w_{1} and w_{2} according to equation .
Pyramidalization angle
The pyramidalization angle is defined via the vectors v_{ba} and
The pyramidalization angle is then given as 90^{°} − θ(v_{ba},w_{1}).
CremerPople parameters
The definitions of the CremerPople parameters for 5 and 6membered rings is described in [40].
Boeyens classification
The Boeyens classification scheme is described in [41].
8.7 Kinetic energy adjustments
There are several options how to adjust the kinetic energy after a surface hop occured. The simplest option is to not adjust the kinetic energy at all (input: ekincorrect none), but this obviously leads to violation of the conservation of total energy.
Alternatively, the velocities of all atoms can be rescaled so that the new kinetic energy and the potential energy of the new state β again sum up to the total energy.
Within this methodology, when the energy of the old state α was lower than the energy of the new state β the kinetic energy is lowered. Since the kinetic energy must be positive, this implies that there might be states which cannot be reached (their potential energy is above the total energy). A hop to such a state is called “frustrated hop” and will be rejected by SHARC. Rescaling parallel to the nuclear velocity vector is requested with ekincorrect parallel_vel.
Alternatively, according to Tully’s original formulation of the surface hopping method [14], after a hop from α to β only the component of the velocity along the direction of the nonadiabatic coupling vector K_{βα} should be rescaled. With
the available energy can be calculated:
If ∆ < 0, the hop is frustrated and will be rejected. Otherwise, the scaled velocities v′ can be calculated as
This procedure can be requested with ekincorrect parallel_nac. Note that in this case SHARC will request the nonadiabatic coupling vectors, even if they are not used for the wavefunction propagation.
8.8 Laser fields
The program laser.x can calculate laser fields as superpositions of several analytical, possibly chirped, laser pulses. In the following, the laser parametrization is given (see [51] for further details).
8.8.1 Form of the laser field
In general, the laser field ϵ(t) is a linear superposition of a number of laser pulses l_{i}(t):
where p_{i} is the normalized polarization vector of pulse i.
A pulse l(t) is formed as the product of an envelope function and a field function.
8.8.2 Envelope functions
There are two types of envelope function defined in laser.x, which are Gaussian and sinusoidal.
The Gaussian envelope is defined as:
where E_{0} is the peak field strength, FWHM is the full width at half maximum and t_{c} is the temporal center of the pulse.
The sinusoidal envelope is defined as:
where again E_{0} is the peak field strength, t_{0} and t_{c} define the interval where the field strength increases, and t_{c2} and t_{e} define the interval where the field strength decreases. Figure 8.1 shows the general form of the envelope functions and the meaning of the temporal parameters t_{0}, t_{c}, t_{2c} and t_{e}.
8.8.3 Field functions
The field function f(t) is defined as:
where ω_{0} is the central frequency and ϕ is the phase of the pulse. Even though the laser field is complex in this expression, in the propagation of the electronic wavefunction in SHARC only the real part is used.
8.8.4 Chirped pulses
In order to apply a chirp to the laser pulse l(t), it is first Fourier transformed to the frequency domain, giving the function ~l(ω). The chirp is applied by calculating:
The chirped laser in the time domain l′(t) is then obtained by Fourier transform of the chirped pulse ~l′(ω).
8.8.5 Quadratic chirp without Fourier transform
If laser.x was compiled without the FFTW package, the only accessible chirps are quadratic chirps for Gaussian pulses:
Other chirps are only possible with the Fourier transformation.
8.9 Laser interactions
The laser field ϵ is included in the propagation of the electronic wavefunction. In each substep of the propagation, the interaction of the laser field with the dipole moments is included in the Hamiltonian. The contribution V_{i} is in each timestep added to the Hamiltonian in equations or , respectively:
where i, n t and ∆t are defined as in section 8.19.
8.9.1 Surface Hopping with laser fields
If laser fields are present, there can be two fundamentally different types of hops: laserinduced hops and nonadiabatic hops. The latter ones are the same hops as in the laserfree simulations, and demand that the total energy is conserved. The laserinduced hops on the other hand demand that the momentum (kinetic energy) is conserved. Hence, SHARC needs to decide for every hop whether it is laserinduced or not.
Consider a previous state α and a new state β. Currently, the hop is classified based on the energy gap ∆E=E_{β}^{diag}−E_{α}^{diag} and the instantaneous central energy of the laser pulse ω.
The hop is assumed to be laserinduced if
where W is a fixed parameter. W can be set using the input keyword laserwidth.
If a hop has been classified as laserfree, the momentum is adjusted according to the equations given in section 8.7.
8.10 Random initial velocities
Random initial velocities are calculated with a given amount of kinetic energy E per atom a. For each atom, the velocity is calculated as follows, with two uniform random numbers θ and ϕ, from the interval [0,1[:
This procedure gives a uniform probability distribution on a sphere with radius √{2E/m_{a}}.
Note that the translational and rotational components of random initial velocities are not projected out in the current implementation.
Random initial velocities can be requested in the input with veloc random E, where E is a float defining the kinetic energy per atom (in eV).
8.11 Representations
Within SHARC, two different representations for the electronic states are used. The first is the socalled MCH basis, which is the basis of the eigenfunctions of the molecular Coulomb Hamiltonian. The molecular Coulomb Hamiltonian is the standard electronic Hamiltonian employed by the majority of quantum chemistry programs. It contains only the kinetic energy of the electrons and the potential energy arising from the Coulomb interaction between the electrons and nuclei.
With this hamiltonian, states of the same multiplicity couple via the nonadiabatic couplings, while states of different multiplicity do not interact at all.
The second representation used in SHARC is the socalled diagonal representation. It is the basis of the eigenfunctions of the total Hamiltonian.
The term ∧H_{el}^{coup} contains additional couplings not contained in the molecular Coulomb Hamiltonian. The most common couplings are spinorbit couplings and interactions with an external electric field.
Both of these couplings introduce offdiagonal elements in the total Hamiltonian. Thus, the eigenfunctions of the molecular Coulomb Hamiltonian are not the eigenfunctions of the total Hamiltonian.
Within SHARC, usually quantum chemistry information is read in the MCH representation, while the surface hopping is performed in the diagonal one.
8.11.1 Current state in MCH representation
Oftentimes, it is very useful to know to which MCH state the currently active diagonal state corresponds. If ∧H_{el}^{coup} is small or the state separation is large, then each diagonal state approximately corresponds to one MCH state. Only in the case of large couplings and/or neardegenerate states are the MCH states strongly mixed in the diagonal states.
In order to obtain for a given timestep from the currently active diagonal state β the corresponding MCH state α, a vector c^{diag} with c_{i}^{diag}=δ_{iβ} is generated. The vector is transformed into the MCH representation
The corresponding MCH state α is the index of the (absolute) largest element of vector c^{MCH}.
8.12 Sampling from Wigner Distribution
The sampling is based on references [37,[38].
The frequency calculation provides a set of vibrational frequencies {ν_{i}} and the corresponding normal mode vectors {n_{i}}, where i runs from 1 to N=3n_{atom}. It also provides the equilibrium geometry R_{eq}.
In order to create an initial condition (R,v), the following procedure is applied. Initially, R_{0}=R_{eq} and v_{0}=0. Then, for each normal mode i, two random numbers P_{i} and Q_{i} are chosen uniformly from the interval [−3,3]. The value of a ground state quantum Wigner distribution for these values is calculated:
W_{i} is compared to a uniform random r_{i} number from [0,1]. If W_{i} > r_{i}, then P_{i} and Q_{i} are accepted and the coordinates and velocities are updated:
The random number procedure and updates are repeated for all normal modes, until (R_{N},v)_{N} is obtained, which constitutes one initial condition. Finally, the center of mass is restored and translational and rotational components are projected out of v. The harmonic potential energy is given by:
8.13 Scaling
The scaling factor (keyword scaling) applies to all energies and derivatives of energies. Hence, the full Hamiltonian is scaled, and the gradients are scaled. Nothing else is scaled (no dipole moments, nonadiabatic couplings, overlaps, etc).
8.14 Seeding of the RNG
The standard Fortran 90 random number generator is seeded by a sequence of integers of length n, where n depends on the computer architecture. The input of SHARC, however, takes only a single RNG seed, which must reproducibly produce the same sequence of random numbers for the same input.
In order to generate the seed sequence from the single input x, the following procedure is applied:
 Query for the number n,
 Generate a first seed sequence s with s_{i}=x+37i+17i^{2},
 Seed with the sequence s,
 Obtain a sequence r of n random numbers on the interval [0,1[,
 Generate a second seed sequence s′ with s_{i}′=int(65536(r_{i}−[1/2])),
 Reseed with the sequence s′.
The fifth step will generate a sequence of nearly uncorrelated numbers, distributed uniformly over the full range of possible integer values.
8.15 Selection of gradients and nonadiabatic couplings
In order to increase performance, it is possible to omit the calculation of certain gradients and nonadiabatic couplings. An energygapbased algorithm selects at each timestep a subset of all possible gradients and nonadiabatic couplings to be calculated. Given the diagonal energy E^{diag}_{ξ} of the current active state ξ, the gradient g^{MCH}_{α} of MCH state α is calculated if:
where ε_{grad} is the selection threshold.
Similarly, a nonadiabatic coupling vector K^{MCH}_{βα} is calculated if:
with selection threshold ε_{nac}.
Neither g^{MCH}_{α} nor K^{MCH}_{βα} are ever calculated if α or β are frozen states.
There is only one keyword (eselect) to set the selection threshold, so ε_{grad} and ε_{nac} are the same in most cases.
8.16 State ordering
The canonical ordering of MCH states of different S and M_{S} in SHARC is as follows. In the innermost loop, the quantum number is increased; then M_{S} and finally S. Example:
nstates 3 0 3
In this example, the order of states is given as:
Number  Label  S  M_{S}  n 
1  S_{0}  0  0  1 
2  S_{1}  0  0  2 
3  S_{2}  0  0  3 
4  T_{1}^{−}  2  1  1 
5  T_{2}^{−}  2  1  2 
6  T_{3}^{−}  2  1  3 
7  T_{1}^{0}  2  0  1 
8  T_{2}^{0}  2  0  2 
9  T_{3}^{0}  2  0  3 
10  T_{1}^{+}  2  +1  1 
11  T_{2}^{+}  2  +1  2 
12  T_{3}^{+}  2  +1  3 
The canonical ordering of states is for example important in order to specify the initial state in the MCH basis (using the state keyword in the input file).
8.17 Surface Hopping
Given two coefficient vectors c^{diag}(t) and c^{diag}(t+∆t) and the corresponding propagator matrix R^{MCH}(t+∆t,t), the surface hopping probabilities are given by
where, however, P_{β→β}=0 and all negative P_{β→α} are set to zero.
The hopping procedure itself obtains a uniform random number r from the interval [0,1]. A hop to state α is performed, if
See section 8.7 for further details on how frustrated hops are handled.
8.18 Velocity Verlet
The nuclear coordinates of atom A are updated according to the Velocity Verlet algorithm [52], based on the gradient of state β at R(t) and R(t+∆t):
Currently, there are no other integrators for the nuclear motion implemented in SHARC.
8.19 Wavefunction propagation
The electronic wavefunction is needed in order to carry out surface hopping. The electronic wavefunction is expanded in the basis of the socalled model space S, which includes the few lowest states ψ^{MCH}_{α}〉 of the multiplicities under consideration (e.g. the 3 lowest singlet and 2 lowest triplet states).
All multiplet components are included explicitly, i.e., the inclusion of an MCH triplet state adds three explicit states to the model space (the three components of the triplet).
Within SHARC, the wavefunction is represented just by the vector c^{MCH}. The Hamiltonian H^{MCH} is represented in matrix form with elements:
From the MCH representation, the diagonal representation can be obtained by unitary transformation within the model space S (U^{†}H^{MCH}U=H^{diag} and U^{†}c^{MCH}=c^{diag}):
and
The propagation of the electronic wavefunction from time t to t+∆t can then be written as the product of a propagation matrix with the coefficients at time t:
or
In order to calculate R^{MCH}(t+∆t,t), SHARC uses (unitary) operator exponentials.
8.19.1 Propagation using nonadiabatic couplings
Here we assume that in the dynamics the interaction between the electronic states is described by a matrix of nonadiabatic couplings K^{MCH}(t), such that
or
In equation , the timederivative couplings are directly calculated by the quantum chemistry program (use coupling ddt in the SHARC input), while in the matrix K^{MCH}(t) is obtained from the scalar product of the nuclear velocity and the nonadiabatic coupling vectors (use coupling ddr in the input).
The propagation matrix can then be written as
with the timeordering operator ∧T. For small timesteps ∆t, H^{MCH}(τ) and K^{MCH}(τ) can be interpolated linearly
And in order to have a sufficiently small timestep for this to work, the interval (t,t+∆t) is further split into subtimesteps ∆τ = [(∆t)/n].
8.19.2 Propagation using overlap matrices
In many situations, the nonadiabatic couplings in K^{MCH} are very localized on the potential hypersurfaces. If this is the case, in the dynamics very short timesteps are necessary to properly sample the nonadiabatic couplings. If too large timesteps are used, part of the coupling may be missed, leading to wrong population transfer. The local diabatization algorithm gives more numerical stability in these situations. It can be requested with the line coupling overlap in the input file.
Within this algorithm, the change of the electronic states between timesteps is described by the overlap matrix S^{MCH}(t,t+∆t)
With this, the propagator matrix can be written as
Chapter 9 Bibliography
Bibliography
 [1]

M. Kasha:
“Characterization of
electronic transitions in complex molecules”. Discuss. Faraday Soc.,
9, 14 (1950).  [2]

C. M. Marian: “Spinorbit
coupling and intersystem crossing in molecules”. WIREs Comput. Mol.
Sci., 2, 187203 (2012).  [3]

I. Wilkinson, A. E. Boguslavskiy, J. Mikosch, D. M. Villeneuve, H.J. Wörner,
M. Spanner, S. Patchkovskii, A. Stolow:
“Nonadiabatic and
intersystem crossing dynamics in SO_{2} I: Bound State Relaxation Studied
By TimeResolved Photoelectron Photoion Coincidence Spectroscopy”. J.
Chem. Phys., 140, 204 301 (2014).  [4]

S. Mai, P. Marquetand, L. González:
“NonAdiabatic Dynamics in
SO_{2}: II. The Role of Triplet States Studied by SurfaceHopping
Simulations”. J. Chem. Phys., 140, 204 302 (2014).  [5]

C. Lévêque, R. Taïeb, H. Köppel:
“Communication:
Theoretical prediction of the importance of the ^{3}B_{2} state in the dynamics
of sulfur dioxide”. J. Chem. Phys., 140, 091101 (2014).  [6]

T. J. Penfold, R. Spesyvtsev, O. M. Kirkby, R. S. Minns, D. S. N. Parker, H. H.
Fielding, G. A. Worth:
“Quantum dynamics study of
the competing ultrafast intersystem crossing and internal conversion in the
“channel 3″ region of benzene”. J. Chem. Phys., 137,
204310 (2012).  [7]

R. A. Vogt, C. Reichardt, C. E. CrespoHernández:
“ExcitedState Dynamics in
NitroNaphthalene Derivatives: Intersystem Crossing to the Triplet Manifold
in Hundreds of Femtoseconds”. J. Phys. Chem., 117,
65806588 (2013).  [8]

C. E. CrespoHernández, B. Cohen, P. M. Hare, B. Kohler:
“Ultrafast ExcitedState
Dynamics in Nucleic Acids”. Chem. Rev., 104, 19772020
(2004).  [9]

M. Richter, P. Marquetand, J. GonzálezVázquez, I. Sola, L. González:
“Femtosecond Intersystem
Crossing in the DNA Nucleobase Cytosine”. J. Phys. Chem. Lett.,
3, 30903095 (2012).  [10]

L. MartínezFernández, L. González, I. Corral:
“An Ab Initio Mechanism
for Efficient Population of Triplet States in Cytotoxic Sulfur Substituted
DNA Bases: The Case of 6Thioguanine”. Chem. Commun., 48,
21342136 (2012).  [11]

S. Mai, P. Marquetand, M. Richter, J. GonzálezVázquez, L. González:
“Singlet and Triplet
ExcitedState Dynamics Study of the Keto and Enol Tautomers of Cytosine”.
ChemPhysChem, 14, 29202931 (2013).  [12]

C. Reichardt, C. E. CrespoHernández:
“Ultrafast Spin Crossover
in 4Thiothymidine in an Ionic Liquid”. Chem. Commun., 46,
59635965 (2010).  [13]

J. C. Tully, R. K. Preston:
“Trajectory Surface
Hopping Approach to Nonadiabatic Molecular Collisions: The Reaction of
H^{+} with D_{2}“. J. Chem. Phys., 55, 562572 (1971).  [14]

J. C. Tully: “Molecular
dynamics with electronic transitions”. J. Chem. Phys., 93,
10611071 (1990).  [15]

M. Barbatti: “Nonadiabatic
dynamics with trajectory surface hopping method”. WIREs Comput. Mol.
Sci., 1, 620633 (2011).  [16]

N. L. Doltsinis: Molecular Dynamics Beyond the BornOppenheimer
Approximation: Mixed QuantumClassical Approaches, volume 31 of NIC
Series. John von Neuman Institut for Computing (2006).  [17]

N. L. Doltsinis, D. Marx:
“First Principles
Molecular Dynamics Involving Excited States And Nonadiabatic Transitions”.
J. Theor. Comput. Chem., 1, 319349 (2002).  [18]

S. HammesSchiffer, J. C. Tully:
“Proton transfer in
solution: Molecular dynamics with quantum transitions”. J. Chem.
Phys., 101, 46574667 (1994).  [19]

G. Granucci, M. Persico:
“Critical appraisal of the
fewest switches algorithm for surface hopping”. J. Chem. Phys.,
126, 134 114 (2007).  [20]

M. Thachuk, M. Y. Ivanov, D. M. Wardlaw:
“A semiclassical approach
to intensefield abovethreshold dissociation in the long wavelength limit”.
J. Chem. Phys., 105, 40944104 (1996).  [21]

B. Maiti, G. C. Schatz, G. Lendvay:
“Importance of Intersystem
Crossing in the S(3P, 1D) + H2 → SH + H Reaction”. J.
Phys. Chem. A, 108, 87728781 (2004).  [22]

G. A. Jones, A. Acocella, F. Zerbetto:
“OntheFly,
ElectricFieldDriven, Coupled ElectronNuclear Dynamics”. J. Phys.
Chem. A, 112, 96509656 (2008).  [23]

R. Mitri\’c, J. Petersen, V. Bonac\’cKoutecký:
“Laserfieldinduced
surfacehopping method for the simulation and control of ultrafast
photodynamics”. Phys. Rev. A, 79, 053 416 (2009).  [24]

G. Granucci, M. Persico, G. Spighi:
“Surface hopping
trajectory simulations with spinorbit and dynamical couplings”. J.
Chem. Phys., 137, 22A501 (2012).  [25]

B. F. E. Curchod, T. J. Penfold, U. Rothlisberger, I. Tavernelli:
“Local Control
Theory using Trajectory Surface Hopping and LinearResponse TimeDependent
Density Functional Theory”. Chimia, 67, 218221 (2013).  [26]

G. Cui, W. Thiel:
“Generalized
trajectory surfacehopping method for internal conversion and intersystem
crossing”. J. Chem. Phys., 141, 124 101 (2014).  [27]

J. Pittner, H. Lischka, M. Barbatti:
“Optimization
of mixed quantumclassical dynamics: Timederivative coupling terms and
selected couplings”. Chem. Phys., 356, 147 – 152 (2009).  [28]

M. Richter, P. Marquetand, J. GonzálezVázquez, I. Sola, L. González:
“SHARC: ab Initio
Molecular Dynamics with Surface Hopping in the Adiabatic Representation
Including Arbitrary Couplings”. J. Chem. Theory Comput., 7,
12531258 (2011).  [29]

S. Mai, M. Richter, M. Ruckenbauer, M. Oppel, P. Marquetand, L. González:
“SHARC: Surface Hopping Including Arbitrary Couplings – Program
Package for NonAdiabatic Dynamics”. sharcmd.org (2014).  [30]

M. Richter, P. Marquetand, J. GonzálezVázquez, I. Sola, L. González:
“Correction to SHARC: Ab
Initio Molecular Dynamics with Surface Hopping in the Adiabatic
Representation Including Arbitrary Couplings?” J. Chem. Theory
Comput., 8, 374374 (2012).  [31]

J. J. Bajo, J. GonzálezVázquez, I. Sola, J. Santamaria, M. Richter,
P. Marquetand, L. González:
“Mixed QuantumClassical
Dynamics in the Adiabatic Representation to simulate molecules driven by
strong laser pulses”. J. Phys. Chem. A, 116, 28002807
(2011).  [32]

P. Marquetand, M. Richter, J. GonzálezVázquez, I. Sola, L. González:
“Nonadiabatic ab initio
molecular dynamics including spinorbit coupling and laser fields”.
Faraday Discuss., 153, 261273 (2011).  [33]

S. Mai, M. Richter, P. Marquetand, L. González:
“Excitation of
Nucleobases from a Computational Perspective II: Dynamics”. In
M. Barbatti, A. C. Borin, S. Ullrich (editors), Photoinduced Phenomena
in Nucleic Acids, Topics in Current Chemistry, Springer Berlin Heidelberg
(2014).  [34]

L. González, P. Marquetand, M. Richter, J. GonzálezVázquez, I. Sola:
“Ultrafast
laserinduced processes described by ab initio molecular dynamics”. In
R. de Nalda, L. Bañares (editors), Ultrafast Phenomena in Molecular
Sciences, volume 107 of Springer Series in Chemical Physics,
145170, Springer (2014).  [35]

G. Granucci, M. Persico, A. Toniolo:
“Direct semiclassical
simulation of photochemical processes with semiempirical wave functions”.
J. Chem. Phys., 114, 10 60810 615 (2001).  [36]

F. Plasser, G. Granucci, J. Pittner, M. Barbatti, M. Persico, H. Lischka:
“Surface hopping
dynamics using a locally diabatic formalism: Charge transfer in the ethylene
dimer cation and excited state dynamics in the 2pyridone dimer”. J.
Chem. Phys., 137, 22A51422A51413 (2012).  [37]

J. P. Dahl, M. Springborg:
“The Morse oscillator
in position space, momentum space, and phase space”. J. Chem. Phys.,
88, 45354547 (1988).  [38]

R. Schinke: Photodissociation Dynamics: Spectroscopy and Fragmentation of
Small Polyatomic Molecules. Cambridge University Press (1995).  [39]

M. Barbatti, G. Granucci, M. Ruckenbauer, F. Plasser, J. Pittner, M. Persico,
H. Lischka: “NEWTONX: a package for Newtonian dynamics close to
the crossing seam, version 1.2”. www.newtonx.org (2011).  [40]

D. Cremer, J. A. Pople:
“General definition of
ring puckering coordinates”. J. Am. Chem. Soc., 97,
13541358 (1975).  [41]

J. C. A. Boeyens: “The
conformation of sixmembered rings”. J. Cryst. Mol. Struct.,
8, 317 (1978).  [42]

H. Lischka, T. Müller, P. G. Szalay, I. Shavitt, R. M. Pitzer, R. Shepard:
“Columbus – a program
system for advanced multireference theory calculations.” WIREs Comput.
Mol. Sci., 1, 191199 (2011).  [43]

H. Lischka, R. Shepard, I. Shavitt, R. M. Pitzer, M. Dallos, T. Müller, P. G.
Szalay, F. B. Brown, R. Ahlrichs, H. J. Böhm, A. Chang, D. C. Comeau,
R. Gdanitz, H. Dachsel, C. Ehrhardt, M. Ernzerhof, P. Höchtl, S. Irle,
G. Kedziora, T. Kovar, V. Parasuk, M. J. M. Pepper, P. Scharf, H. Schiffer,
M. Schindler, M. Schüler, M. Seth, E. A. Stahlberg, J.G. Zhao,
S. Yabushita, Z. Zhang, M. Barbatti, S. Matsika, M. Schuurmann, D. R.
Yarkony, S. R. Brozell, E. V. Beck, , J.P. Blaudeau, M. Ruckenbauer,
B. Sellner, F. Plasser, J. J. Szymczak: “COLUMBUS, an ab initio
electronic structure program, release 7.0” (2012),
http://www.univie.ac.at/columbus/docs_COL70/documentation_main.html.  [44]

S. Yabushita, Z. Zhang, R. M. Pitzer:
“SpinOrbit
Configuration Interaction Using the Graphical Unitary Group Approach and
Relativistic Core Potential and SpinOrbit Operators”. J. Phys.
Chem. A, 103, 57915800 (1999).  [45]

S. Mai, T. Müller, P. Marquetand, F. Plasser, H. Lischka, L. González:
“Perturbational Treatment
of SpinOrbit Coupling for Generally Applicable HighLevel MultiReference
Methods”. J. Chem. Phys., 141, 074 105 (2014).  [46]

H. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schütz:
“Molpro: a generalpurpose
quantum chemistry program package”. WIREs Comput. Mol. Sci.,
2, 242253 (2012).  [47]

H.J. Werner, P. J. Knowles, G. Knizia, F. R. Manby, M. Schütz, et al.:
“MOLPRO, version 2012.1, a package of ab initio programs”.
www.molpro.net (2012).  [48]

G. Karlström, R. Lindh, P.Å. Malmqvist, B. O. Roos, U. Ryde,
V. Veryazov, P.O. Widmark, M. Cossi, B. Schimmelpfennig, P. Neogrady,
L. Seijo:
“MOLCAS: a
program package for computational chemistry”. Comput. Mater. Sci.,
28, 222 – 239 (2003), proceedings of the Symposium on Software
Development for Process and Materials Design.  [49]

F. Aquilante, L. De Vico, N. Ferré, G. Ghigo, P.Å. Malmqvist,
P. Neogrády, T. B. Pedersen, M. Pitonák, M. Reiher, B. O. Roos,
L. SerranoAndrés, M. Urban, V. Veryazov, R. Lindh:
“MOLCAS 7: The Next
Generation”. J. Comput. Chem., 31, 224247 (2010),
www.molcas.org.  [50]

G. Granucci, M. Persico, A. Zoccante:
“Including quantum
decoherence in surface hopping”. J. Chem. Phys., 133,
134 111 (2010).  [51]

P. Marquetand: Vectorial properties and laser control of molecular
dynamics. Ph.D. thesis, University of Würzburg (2007),
http://opus.bibliothek.uniwuerzburg.de/volltexte/2007/2469/.  [52]

L. Verlet: “Computer
Ëxperiments” on Classical Fluids. I. Thermodynamical Properties of
LennardJones Molecules”. Phys. Rev., 159, 98103 (1967).