QCMagic

QCMagic is a set of Python libraries and scripts to manipulate Q-Chem output files. To get access to QCMagic, visit

   https://gitlab.developers.cam.ac.uk/ch/thom/qcmagic

and log in with your Raven account. If you do not have permissions to this site, please send an email to Alex Thom to request access. Once you have logged on to the GitLab Repository, go to Settings -> SSH Keys and follow the instructions to add your own public key to the Repository. Once done, you will be able to clone QCMagic with

   git clone git@gitlab.developers.cam.ac.uk:ch/thom/qcmagic.git

The current main branch is master. If you wish to develop QCMagic, you need to create your own branch using

   git checkout -b yourNewBranch

and then verify that you are indeed on your own branch:

   git status

Usage

Having cloned QCMagic, assuming it is in your home directory, add the following to your .bashrc:

  export PATH=$PATH:$HOME/qcmagic
  export PYTHONPATH=$PYTHONPATH:$HOME/qcmagic
  export VMDSIMPLEGRAPHICS=1

The last is to ensure you can open a VMD window (provided X forwarding is enabled, and you have VMD installed locally).

If there are problems using VMD via SSH try using xpra

Probably the best way to explore the capabilities of qcmagic is to use iPython. Say you have a Q-Chem output file myfile.out containing a single job; in iPython type:

  from QCManager import *
  O = OutputFile('myfile.out') # The class OutputFile is defined in QCOutput.py
  O.Parse()                    # Extract data from myfile.out. Returns False if no Q-Chem errors are found.

Some data will be stored as values in a container, which is a special type of dictionary defined in QCSupport.py, whilst others are regular attributes of O. For example, if myfile.out is the result of a Q-Chem optimisation ('jobtype opt'), you could find the initial and final energies with:

  Ei = O.Steps[0]['InternalEnergy']
  Ef = O.Steps[-1]['InternalEnergy']

and the initial and final geometries with:

  Xi = O.Steps[0]['System'].Atoms
  Xf = O.Steps[-1]['System'].Atoms
  # Or equivalently:
  Xf = O.System.Atoms

makePlots.py

Note This currently appears to work with qchem input files, but not output files.

Example command:

  python makePlots.py my_qchem_file.inp -T tag

This runs a qchem calculation on my_qchem_file.inp and plots the orbitals in VMD. Output files will be saved to a new directory tag, or Unnamed if tag is unspecified. Run with --help to view all possible options.

Once the plot is up, the display can be controlled entirely from the terminal. For example, view orbital #3 by typing 3. Type help to view all available commands.

scanSurface.py

Overview

A general command for scanSurface has the following form:

  scanSurface.py [options] qchem_file [qchem_file_2 ...] Tag

where qchem_file's are old Q-Chem input or output files whose inputs are used as the system (and assumed to be setup as above). If these files contain saved minima, they are stored and relocated. Outputs from scanSurface are stored in Tag.sd (a pickled Python data structure containing all calculation results), Tag.dat (a printout of energies at different geometries), and when necessary, Tag.dat2 (a printout of NOCI energies at different geometries).

The principle of working of scanSurface is rather simple: from the initial geometry of the system supplied in the old Q-Chem inputs or outputs, each electronic state requested will be followed as the geometry of the molecule is varied. The electron density of the state at the previous geometry is read in as the initial guess for the convergence procedure performed by Q-Chem at the current geometry. Upon convergence, a state distance is calculated between the electron densities at the current and previous geometry, so that only when the new density is located within a certain maximum threshold of the old density can it be accepted as belonging to the same state. scanSurface saves the Q-Chem output file for every calculation as Run*.out.

An example command is:

  scanSurface.py --stretch=0,1,1,0.01 my_qchem_file.out Tag

which runs a scan on a molecule specified in the input section of my_qchem_file.out, calling Q-Chem at points on a surface defined by the stretching of the bond between atom0 and atom1 by 10 units in steps of 0.01, and writing the results to Tag.sd and Tag.dat.

If multiple options are used, such as

  scanSurface.py --stretch=0,1,1,0.01 --stretch=0,2,1,0.01 my_qchem_file.out Tag

the scan will zigzag to cover the whole surface.

Options

scanSurface can be run with --help to view all possible options. In this Section, important remarks for some options are presented.

--read-minima

This option is used to specify the minima located by metadynamics in the specified Q-Chem output file that scanSurface should read in and follow. It must be noted that this is one-based and not zero-based, so --read-minima=1 ensures only the first minimum located by metadynamics is scanned. On the other hand, --read-minima=0 would have no effect on the scan. For example,

  scanSurface.py --read-minima=1,2,6 --stretch=0,1,1,0.01 my_qchem_file.out Tag

instructs scanSurface to read in minima 1, 2 and 6 from the old Q-Chem output file my_qchem_file.out, follow each of them as the bond between atom0 and atom1 is stretched by 10 units in steps of 0.01, then store the results as states 0, 1 and 2 in Tag.sd and print out the energies in Tag.dat.

Note that in order for scanSurface to read in the requested minima, the Q-Chem output file must have been generated with the rem option print_orbitals set to true.

--read-only

If on, scanSurface only reads and stores the states from the Q-Chem output file to an SD file without trying to reconverge to any of them. This is useful when one seeks to export the states located by metadynamics in a Q-Chem output file to an SD file for further manipulation, analysis or extraction.

--normal-rotate

This option may be used in simple cases to change a single bond angle. It relies on being able to work out the connectivity from knowledge of the bond lengths and VDW radii, so can be problematic if (a) the bond lengths are unusual, and (b) there are multiple paths between two atoms. For example,

   --normal-rotate=0,1,2,90,0.1

will:

• split the molecule between atom0 and atom1;
• check if each of the other atoms is connected to atom0; if they are, group them with atom0 and if not, group them with atom1; and
• rotate the group of atoms containing atom1 about an axis through atom1, and normal to the plane containing atom0, atom1 and atom2.

--tie

This option allows multiple parameters to change in sync with one another. For example, one might want to investigate the effect of varying two bond lengths at the same time, but would like to ignore what happens when one is stretched and the other unchanged. This could be very useful in the context of molecular vibrations, or when certain symmetry elements need to be preserved.

For each parameter that needs to be tied, two integers are required:

• one referring to the option:

ext_charge -- 0, scale -- 1, stretch -- 2, rotate -- 3, normal_rotate -- 4, orth_rotate -- 5, all_normal_rotate -- 6, rotate_axis -- 7;

• another one specifying the particular instance of said option (zero-based).

An example of a command using --tie:

  scanSurface.py --stretch=0,1,1,0.01 --rotate=2,3,90,0.9 --scale=2 --tie=2,0,3,0 my_qchem_file.out Tag

This will:

• simultaneously stretch bond atom0-atom1 by 1 unit, and rotate about bond atom2-atom3 by 90°, in 100 steps;
• scale the whole molecule up to twice its original size. As the --scale option is not tied, scanSurface will scale separately every possible geometry allowed by the previous (tied) options.

Note that an equal number of steps for each tied parameter must be ensured in order for --tie to work correctly!

--focus

If this flag is on, once scanSurface fails to locate the correct minimum at a geometry, it will not attempt to find other minima at that geometry but move on to the next geometry instead, where it will use the minimum on the same state at the last well-converged geometry as reference. If scanSurface locates the correct minimum after a few other minima that belong to different states, only the correct minimum is saved and used as reference for the next geometry, while all other superfluous minima are discarded.

This option is useful when following states that pass through regions of discontinuity (and non-convergence), as it prevents scanSurface from getting trapped in a cascading landslide where it tries to look for other minima and follow them, only to be stuck in more regions of discontinuity (and non-convergence) of these new minima later on, halting the state-tracking process altogether.

qcSDExtract.py

Currently qcSDExtract only works with traditional, non-augmented SDS objects (i.e. those that do not contain variations in exchange-correlation).

Overview

qcSDExtract reads in states stored in a single SD file generated by scanSurface, performs further calculations on them if need be, and extracts out their energies, atomic charges and spins from Mulliken analysis, as well as oxidation states of individual atoms from LOBA analysis. qcSDExtract can parametrise geometries using multiple bond lengths, bond angles, and dihedral angles, allowing multi-dimensional surfaces to be traced out.

A general command for qcSDExtract has the following form:

  qcSDExtract.py [options] SDfile Tag

where SD_file is an SD file containing information of all states read in and possibly followed by scanSurface. Several output files are generated by qcSDExtract:

• Tag.energy, containing all geometrical parameters requested and the energies of the associated states;
• Tag.atomi.charge, containing all geometrical parameters requested and the Mullikan charges of atomi in the associated states;
• Tag.atomi.spinz, containing all geometrical parameters requested and the Mullikan spin charges (i.e. spin imbalances) of atomi in the associated states; and
• Tag.atomi.oxstate, containing all geometrical parameters requested and the LOBA oxidation states of atomi in the associated states.

An example command is:

  qcSDExtract.py --read-minima=0,1,2 --read-geometries=0,1,2,3,4 --bondlength=0,1,0,2,1,2 --bondangle=1,0,2 --rem="LOBA 12" --atoms=0 mydata.sd Tag

which:

• reads in states 0, 1 and 2 at geometries 0, 1, 2, 3 and 4 stored in my_data.sd;
• parametrises these geometries using bond lengths atom0--atom1, atom0--atom2, atom1--atom2, and bond angle atom1--atom0--atom2;
• performs additional LOBA 12 calculations on these minima; and
• outputs Tag.energy, Tag.atom0.charge, Tag.atom0.spinz and Tag.atom0.oxstate, where the three latter files contain details for atom0.

In addition, the actual Q-Chem output file generated when qcSDExtract calls Q-Chem to reconverge into state j at geometry i (i and j are 0-based indices) and/or to perform additional calculations is saved as _SDfile.gi.sj.out (note the underscore prefix).

Options

--read-minima

At least one 0-based index of one or more states to read in from the SD file must be specified. Multiple indices are separated by commas. For example,

   --read-minima=0,1,2

instructs qcSDExtract to read in minima 0, 1 and 2 from the SD file.

Q: How do I know which (1-based) minimum in the original Q-Chem output file corresponds to which (0-based) minimum in the SD file after scanSurface has been run?

A: Examine the stdout of scanSurface and look for the pattern "Old min XYZ -> New sol ABC" which shows the mapping of the (1-based) indices in the original Q-Chem output file, XYZ, to the (0-based) indices in the SD file, ABC.

--read-geometries

This option allows one to specify the 0-based indices of one or more geometries to be examined, separated by commas. If no indices are specified, or if this option is omitted altogether, all geometries will be scanned. For example, if only geometries with indices 0, 5 and 10 are needed, one should specify

   --read-geometries=0,5,10

which is useful when only certain interesting portions of the surface are required to be examined.

--atoms

This option allows the specification of a list of atom indices (0-based) from which chemical properties are extracted. If this option is omitted, properties will be extracted from only atom0. For example, to extract chemical properties of atom0 and atom1 along state 0 at all geometries, one would use the command

   qcSDExtract.py --read-minima=0 --atoms=0,1 --rem="LOBA 12" my_file.sd Tag

--bondlength

This is where one can specify the bond lengths that should be used to parametrise the geometries. This, together with the following two options, --bondangle and --bonddihedral, provides a workaround to construct a full multi-dimensional surface from the flatten-down one-dimensional storage of geometries originally implemented in scanSurface.

Here, one specifies a list of 0-based indices of pairs of atoms whose bond lengths will be used to parametrise the geometries. Therefore, the list must contain an even number of elements, otherwise qcSDExtract will exit with an error. For example,

   --bondlength=0,1,0,2,1,2

parametrises the geometries using three bond lengths: atom0--atom1, atom0--atom2, and atom1--atom2. If this option is omitted, atom0--atom1 bond length will be used by default.

--bondangle

Here, one specifies a list of 0-based indices of triplets of atoms whose bond angles will be used to parametrise the geometries. The list must contain a multiple of three of elements. In each triplet, the middle index is the centre atom. For example,

   --bondangle=1,0,2,1,0,3

uses the bond angles atom1--atom0--atom2 and atom1--atom0--atom3 to parametrise the geometries. There is no default for this option.

--bonddihedral

This option specifies a list of 0-based indices of triplets of atoms whose bond dihedrals will be used to parametrise the geometries. The list must contain a multiple of four of elements. In each quartet, the first three atoms form the initial plane. For example,

   --bonddihedral=0,1,2,3

uses the dihedral angle between the two planes atom0--atom1--atom2 and atom1--atom2--atom3 to parametrise the geometries. There is no default for this option.

varyXCB.py

Overview

varyXCB reads in minima from a Q-Chem output file or an SD file (traditional or augmented) and follows them as exchange-correlation functionals are varied, or as basis sets are changed.

A general command for varyXCB has the following form:

  varyXCB.py [options] file Tag

where file is a Q-Chem output file or an SD file containing minima to be read in. If the SD file contains an instance of class SystemDataSet, it will be augmented to an instance of class SystemDataSetXC, placing all existing minima at exchange-correlation functional index 0 (see later). However, if the SD file already contains a SystemDataSetXC object, it will be read in as-is, but an index for exchange-correlation functionals must be specified via --read-xc so that varyXCB will know which minima to follow (see later).

When run in exchange-correlation following mode, varyXCB outputs the energies of the minima along the parametrised path in the exchange-correlation space in Tag.gi.sj.dat, and stores all minima at all exchange-correlation functionals along the path in Tag.xc.sd, which contains a SystemDataSetXC object. On the other hand, when run in basis projection mode, varyXCB outputs the energies of the minima in the new basis set in Tag.newbasis.gi.sj.dat, and stores all minima in the new basis set in Tag.b.sd, which contains a SystemDataSetXC object.

An example command for exchange-correlation following mode is:

  varyXCB.py --xc-template=xctemplate --read-minima=1,2,3 --para=lambda,0,1,0.1 my_qchem_file.out Tag

which:

• reads in minima 1, 2 and 3 from my_qchem_file.out and stores them at exchange-correlation functional index 0 of states 0, 1 and 2 respectively in the SystemDataSetXC object in Tag.xc.sd;
• uses the content of xctemplate as a template for $xc_functional section in Q-Chem runs;
• understands that lambda in xctemplate represents a variable parameter that needs to be increased from 0 to 1 in steps of 0.1 to define a path in the exchange-correlation space;
• runs Q-Chem at different exchange-correlation functionals; and
• stores all calculation results in Tag.xc.sd and outputs Tag.gi.sj.dat for j=0,1,2 which contains the energies of states 1, 2 and 3 (read from the Q-Chem output file) as followed along the specified exchange-correlation path.

Another example command for basis projection mode is:

  varyXCB.py --basis-template=changebasis.template --purecart=1112 --read-minima=0,2,5 --read-geometries=0 --read-xc=1
             --rem="SCF_MAX_CYCLES 100" --rem="SCF_MINFIND_READDISTTHRESH 00500" file.sd Tag

which:

• reads in minima 0, 2 and 5 at geometry 0 and exchange-correlation functional index 1 from file.sd and stores them at exchange-correlation functional index 0 of states 0, 1 and 2 respectively in the SystemDataSetXC object in Tag.b.sd;
• reads in the new basis set specified in changebasis.template;
• uses purecart specification 1112 for the new basis set;
• runs basis projection calculations with additional rem options;
• outputs Tag.newbasis.g0.si.dat containing the energies of state i=0,1,2 at geometry 0 in the old and new basis sets;
• outputs Run-new_basis.gi.sj.xc1.out which is the Q-Chem output file for the actual basis-set projection calculation of state i=0,1,2 at geometry 0 and exchange-correlation 1; and
• stores all calculated data in Tag.b.sd.

Structure of a QCSystemData.SystemDataSetXC Object (for developers only)

QCSystemData.SystemDataSetXC (SystemDataSetXC for short) inherits from the base class QCSystemData.SystemDataSet (SystemDataSet for short) with some existing methods overridden and some new methods implemented to cope with the incorporation of exchange-correlation functional variations.

An SD file is a pickled Python dictionary containing three keys: 'Output', 'Job' and 'SysData'. Let us call this dictionary D. D['SysData'] is then an instance of the class SystemDataSet or the class SystemDataSetXC, depending on how it has been constructed. Here we will only focus on the latter case.

D['SysData'] contains various useful methods and variables. However, all important calculation information is contained within the two variables GeomList and StateList.

GeomList

GeomList is a one-dimensional list containing elements which are instances of class QCMolecule.System. Each element therefore carries attributes for describing and manipulating a particular geometry (i.e. atom arrangement) of the system. For example, GeomList could contain successive geometries of a diatomic molecule, XY, where the bond length X--Y is gradually increased in a particular range.

In other words, if GL is a GeomList, then GL[i] is an instance of class QCMolecule.System and contains information of geometry with index i.

StateList

StateList is a one-dimensional list containing elements which are dictionaries. Each dictionary then contains key-value pairs, each of which corresponds to a state whose the index is given by the key. The corresponding value is yet another list whose elements are instances of QCSystemData.SCFStateRich (a class inheriting from QCSystemData.SCFState and able to store more information such as LOBA oxidation states, atomic charges and spins, as well as exchange-correlation functionals) and describe the state at different exchange-correlation functionals.

In other words, if SL is a StateList, then

• SL is a list with a one-to-one correspondence to GL;
• SL[i] is a dictionary containing information of all states at geometry with index i, GL[i];
• SL[i][j] is a list containing information at all exchange-correlation functionals of state with index j at geometry with index i; and
• SL[i][j][k] is an instance of QCSystemData.SCFStateRich containing information at exchange-correlation functional with index k of state with index j at geometry with index i.

Options

varyXCB can be run with -h or --help to display all options supported and how to use them. Here we explain the two most important options that control the working of varyXCB.

--xc-template

In this option, the path to the template file containing specifications of exchange and correlation functionals in terms of parameters is specified. There is no default value for this option, so if it is omitted, no exchange-correlation tracking will be carried out. The template file must have the following format:

   $xc_functional
       X   exchange_functional_name       parametrised_expression
       ...
       C   correlation_functional_name    parametrised_expression
       ...
       K                                  parametrised_expression !This is an alternative specification for Hartree-Fock exchange -- see Q-Chem manual
   $end

The file must begin with $xc_functional and end with $end, and its format must follow that of the $xc_functional section in Q-Chem inputs, noting that the parametrised expressions must contain no white spaces, e.g.

   1-(0.80*lambda)+gamma

where lambda and gamma are parameters. Any parameters used in the template file must then be declared using the --para option of varyXCB in order for varyXCB to parse the template file correctly. Therefore, for a template file containing the above expression, the following options must be present:

   --para=gamma,0,1,0.5 --para=lambda,0,1,0.1

--basis-template

In this option, the path to the template file containing the specification of a general, larger basis set is indicated. There is no default value for this option, so if it is omitted, no basis set projection will be carried out. This file must begin with $basis and end with $end, and its format must follow that of the $basis section in Q-Chem inputs. In particular,

   $basis
       general_basis_set_specification
   $end

When run in basis projection mode, varyXCB reads in minima in the smaller basis set (this must be one of the standard basis sets in Q-Chem and not a general basis set) and projects them into the larger, general basis set specified in the template file.

Note that the option --purecart MUST be specified when --basis-template is used, so that Q-Chem will know what types of functions (pure or Cartesian) to use for the general basis set.

SymmetryTools

SymmetryTools is a package within QCMagic that allows users to determine the point group of a molecule and use the symmetry elements of the group to generate and analyse symmetry-related and necessarily degenerate single determinants from those already obtained through SCF calculations. SymmetryTools supports proper and improper rotations involving both spatial and spin degrees of freedom. This allows us to generate a complete set of degenerate symmetry-related SCF single determinants and ultimately identify the irreducible representations spanned by them. The full workflow from obtaining multiple SCF solutions using metadynamics to determining the irreducible representations spanned by any one solution and its symmetry partners is shown below.

Step 1: Generating multiple SCF solutions using metadynamics

We begin with obtaining several (say, $N$) SCF solutions using metadynamics in Q-Chem. This involves running an SCF job with SCF_SAVEMINIMA N as a $rem variable where $N\geq 1$. Note that metadynamics is technically when $N>1$. However, when $N=1$, only a single SCF determinant will be saved, just like a normal SCF calculation, but the $rem variable SCF_SAVEMINIMA is needed so that this SCF determinant is outputted in a format that is readable by QCMagic. A sample input named BH3.inp for metadynamics on Q-Chem which prints out 5 energy minima for a borane molecule is shown below.

$molecule
0 1
  B   -3.7312212    1.6982695    0.0000000
  H   -4.8515376    1.3242015    0.0000000
  H   -3.4950153    2.8555260    0.0000000
  H   -2.8471106    0.9150809    0.0000000
$end
$rem
    BASIS STO-3G
    EXCHANGE hf
    SCF_GUESS sad
    SCF_CONVERGENCE 7
    UNRESTRICTED true
    SCF_SAVEMINIMA 5
    SCF_MAX_CYCLES 1000
    MOM_START 1
    PRINT_ORBITALS true
    SCF_MINFIND_INITNORM 09000
    SCF_MINFIND_INITLAMBDA 00200
    SCF_MINFIND_INCREASEFACTOR 10100
    SCF_MINFIND_RANDOMMIXING -15708
    SCF_MINFIND_NRANDOMMIXES 1
    SCF_MINFIND_MIXMETHOD 1
    SCF_MINFIND_MIXENERGY 00100
    SYMMETRY off
    DIIS_SEPARATE_ERRVEC true
    SYM_IGNORE true
$end

We then run

   $ qchem BH3.inp BH3.out

The above job finishes successfully, printing 5 minima in the output file BH3.out. We can briefly examine these minima by their energies. Executing

   $ grep "Saving Minimum" BH3.out

gives

    7     -26.0698458670      3.16E-09  00000 Saving Minimum 1
   60     -23.8075357551      8.98E-08  00000 Saving Minimum 2
  102     -23.5547974285      4.13E-08  00000 Saving Minimum 3
  129     -23.5596201123      9.30E-08  00000 Saving Minimum 4
  145     -25.3032615669      1.80E-08  00000 Saving Minimum 5

Your metadynamics run might find different solutions due to the random nature of metadynamics. Nevertheless, it should still be able to locate five minima fairly easily. We will proceed to determine the irreducible representations spanned by each of these minima together with its symmetry partners. If you wish to separate the minima into different directories, follow step 2 below. If not, you can jump directly to step 3.

Step 2a (optional): Extracting SCF minima from the output file

It might be convenient to extract the minima from the original Q-Chem output file into separate directories for easier management and manipulation. This, however, is not essential for the use of SymmetryTools. To this end, we first need to ask QCMagic to read in the minima from the output file and store them in a .sd file that can be easily re-read in by QCMagic later on. This is done via

   $ runscanSurface.py --read-minima=1,2,3,4,5 --read-only BH3.out BH3-extracted > terminal.out

where

• --read-minima=1,2,3,4,5 instructs runscanSurface.py to read in minima 1--5 from the supplied Q-Chem output file. You can include only the minimum or minima that you are interested in here, and you need not use all minima in the output file. To determine which minima are of interest, you can examine their energies and select those within a certain energy window that you are interested in;
• --read-only instructs runscanSurface.py to only read and store the minima without trying to reconverge into them;
• BH3.out tells runscanSurface.py which output file to read in;
• BH3-extracted specifies the tag name for output files generated by runscanSurface.py (No file extension is needed);
• > terminal.out redirects stdout from runscanSurface.py to terminal.out. This is needed to determine which minima in the .sd file the original minima in the output file have been mapped to.

The above command produces BH3-extracted.dat, BH3-extracted.dat2, BH3-extracted.sd, and terminal.out. BH3-extracted.dat and BH3-extracted.dat2 are empty and need not concern us further since we did not run any reconvergence calculations. BH3-extracted.sd contains the minima extracted from BH3.out. terminal.out contains print statements from runscanSurface.py that are useful for debugging purposes, but most importantly, it contains statements of the form

   Old min $x$ -> New sol $y$

which tells us that minimum $x$ in the Q-Chem output file is now minimum $y$ in the .sd file produced. For example, minima 1--5 in BH3.out have changed to minima 0--4 in BH3-extracted.sd. We then extracts these statements into a file, say oldminnewsol, to be used by a script called runrearrangeSolutions.py later on.

   $ grep "Old min" terminal.out > oldminnewsol

Step 2b (optional): Reconverging into each solution separately

The .sd file generated in step 2a contains all extracted minima. We could in principle use this .sd file for symmetry analysis, but we could also reconverge into each minimum separately, with possible changes to reconvergence criteria (e.g., higher SCF_CONVERGENCE).

As an example, using BH3-extracted.sd from step 2a, we run

   $ runqcSDExtract.py --reconverge --rem="SCF_CONVERGENCE 10" BH3-extracted.sd BH3-highconv

which tries to reconverge into every minimum in BH3-extracted.sd in separate output files. Hence, the number of output files generated is dependent on the number of minima present in the .sd file. In our example, we will have 5 such files:

   BH3-highconv.g0.s$y$.xc0.reconverged.out

where $y$ specifies the solution index number and should corresponds to the $y$ in Old min $x$ -> New sol $y$ obtained earlier.

In addition, the following files are also produced:

• BH3-highconv.xc0.atom0.charge.dat provides the Mulliken net atomic charge on atom 0 for all minima (by default of runqcSDExtract.py);
• BH3-highconv.xc0.atom0.oxstate.dat provides the oxidation state of atom 0 from a LOBA calculation (LOBA is off by default);
• BH3-highconv.xc0.atom0.spinz.dat provides the Mulliken net spinZ ($S_{z}$) on atom 0 (by default of runqcSDExtract.py);
• BH3-highconv.xc0.energy.dat provides the energy of all minima.

Most of the produced .dat files contain information extracted from the various output files and can be used for inspection purposes. At the end of the day, however, only the .out files BH3-highconv.g0.s$y$.xc0.reconverged.out are needed to fully describe the minima as they contain the full set of orbital coefficients.

Step 2c (optional): Separating output files into different directories

It is good practice to analyse each minimum in a separate, properly labelled directory. The script runrearrangeSolutions.py can do this automatically, provided that the files to be moved into a directory called o$x$ contain the string .s$y$. in their names, and that the mapping $x$ -> $y$ is supplied as a line with the format Old min $x$ -> New sol $y$ in the input file passed in to runrearrangeSolutions.py.

For our example, run

   $ runrearrangeSolutions.py oldminnewsol

to arrange each of the BH3-highconv.g0.s$y$.xc0.reconverged.out files into separate directories called o$x$, where the mappings $x$ -> $y$ are provided in oldminnewsol in the correct format specified above.

Step 3: Determining the point group and identifying all symmetry elements of the system

In order to apply symmetry operations on the minima and find the irreducible representations spanned by them, we first need to determine the underlying spatial point group set up by the nuclear framework of the system. This then allows us to identify the symmetry elements present. In our earlier example, since all five solutions come from the same system (i.e., same geometry and basis set), we only need to obtain the point group and symmetry elements from one of them. To this end, run

   $ cd o1
   $ rungeneralisedTransformQCMinima.py foobar BH3-highconv.g0.s0.xc0.reconverged.out BH3-sym

where

• foobar is just a placeholder to satisfy a compulsory positional argument that we do not need or have at the moment;
• BH3-highconv.g0.s0.xc0.reconverged.out is the Q-Chem output file containing all system information, as well as information for state o1;
• BH3-sym specifies the tag name for the files produced by rungeneralisedTransformQCMinima.py.

This produces the following:

• BH3-sym.sym contains an object that has all symmetry and basis information for this system. We will need this object for all symmetry manipulation on all states of this system;
• BH3-sym.g0.sym.summary provides a human-readable summary of the point group and symmetry elements located for this system;
• BH3-sym.sym.sd and BH3-sym.all.transformed.sd do not contain any other useful information for us at the moment.

It is instructive to inspect BH3-sym.g0.sym.summary and check if rungeneralisedTransformQCMinima.py has found a sensible point group and located all expected symmetry elements for the system. In our example, part of this file looks like

*********************************************************************************************************
SYMMETRY INFORMATION
*********************************************************************************************************
Rotational symmetry: Symmetric Top - Oblate, Planar
Point group: D3h
    Distance threshold: 1.00E-02
*********************************************************************************************************  

*********************************************************************************************************
SYMMETRY ELEMENTS LOCATED
*********************************************************************************************************
Proper rotation axes (Cn axes)
Order: 2
Index            x           y           z    alpha/pi     beta/pi    gamma/pi      n,i,power
    0   -0.9485234  -0.3167072   0.0000000   1.0000000   1.0000000  -1.0000000         2,0,p
    1   -0.7485381   0.6630918   0.0000000   0.0000000   1.0000000  -0.0000000         2,1,p
    2   -0.1999852  -0.9797989   0.0000000   0.9359106   1.0000000  -0.9359106         2,2,p 

Order: 3
Index            x           y           z    alpha/pi     beta/pi    gamma/pi      n,i,power
    0    0.0000000   0.0000000  -1.0000000   0.0000000   0.0000000   0.6666667         3,0,p 


Mirror planes (S1 axes)
Index            x           y           z      n,i,power
    0    0.0000000   0.0000000  -1.0000000         1,0,p
    1    0.3167072  -0.9485234   0.0000000         1,1,p
    2   -0.6630918  -0.7485381  -0.0000000         1,2,p
    3    0.9797989  -0.1999852   0.0000000         1,3,p 


Improper rotation axes (Sn axes)
Order: 3
Index            x           y           z      n,i,power
    0    0.0000000   0.0000000  -1.0000000         3,0,p
***********************************************************

We note the following:

• The point group of the system is found to be $D_{3h}$, which is correct.
• Apart from the identity $E$ which is not shown, there are three $C_{2}$ axes, one $C_{3}$ axis, four $\sigma$ planes (one $\sigma _{h}$ and three $\sigma _{v}$), and one $S_{3}$ axis, as expected for $D_{3h}$.
• Each symmetry element is specified by whether it is proper or improper ($C$ or $S$), its order $n$, and its index $i$. These are summarised in the column n,i,power for each symmetry element. The role of the power $p$ will become clear later.
• All symmetry elements are defined with respect to the recentred coordinates where the centre of mass of the molecule has been translated to $(0,0,0)$. Each axis of rotation (proper or improper) passes through the centre of mass and is parallel to the vector whose components are given in the x, y, and z columns.
• Additionally, the Euler angles for the proper rotations are given in $\pi \;\mathrm {rad}$ following the Whitaker convention in the alpha/pi, beta/pi, and gamma/pi columns.

Step 4: Constructing the symmetry operation instruction file

Given a minimum that is a single determinant, we seek to generate other symmetry-related single determinants. To this end, we must first tell rungeneralisedTransformQCMinima.py what symmetry operations it will need to use to generate new single determinants. This is done by means of a file that can contain several lines, each of which contains instructions for several symmetry operations to be applied consecutively from the right to generate a single symmetry-equivalent determinant. More detailed information can be found by running

   $ rungeneralisedTransformQCMinima.py

without any arguments.

For our BH₃ example, however, we shall use a very simple instruction file that contains all spatial symmetry operations (except the trivial identity) in $D_{3h}$ to illustrate how this file should be constructed. Prior to constructing this file, it is important to consult the BH3-sym.g0.sym.summary file generated in step 3 in conjunction with the Character Table for $D_{3h}$. Assuming that that has been done, let us call this file transform.instructions.full whose content is shown below.

# Proper rotations
(C,2,0,1)
(C,2,1,1)
(C,2,2,1)
(C,3,0,1)
(C,-3,0,1)
# Improper rotations
(S,1,0,1) # sigmah
(S,1,1,1) # sigmav
(S,1,2,1) # sigmav
(S,1,3,1) # sigmav
(S,3,0,1)
(S,-3,0,1)

We note the following:

• There are 11 non-comment lines in the file corresponding to the 11 non-identity symmetry operations in $D_{3h}$.
• Anything following # until eol is not considered.
• Each spatial symmetry operation is specified by (signature,order,index,power), where
  • signature gives the type of symmetry operation to be carried out. For our example, signature is either C or S for proper or improper rotations, respectively;
  • order gives the order of rotation. This is given by $n$ in BH3-sym.g0.sym.summary. If $-n$ is used, the inverse of the corresponding rotation with order $n$ is carried out instead;
  • index gives the index of the rotation axis (to distinguish between various axes of the same order. This is given by $i$ in BH3-sym.g0.sym.summary;
  • power gives the number of times the operation is applied.
• Each $C_{2}$ axis gives rise to only one $C_{2}$ operation because, for spatial transformations,$C_{2}=C_{2}^{-1}$. The same applies to mirror planes ($S_{1}$ axes) and inversion centres ($S_{2}$ axes).
• Multiple symmetry operations can be strung together using . and the rightmost symmetry operation is applied first. For example, since $C_{3}^{-1}=C_{2,1}\cdot C_{2,2}$ in our example (try inspecting the effects of these operations based on their axes defined in BH3-sym.g0.sym.summary), we could replace (C,-3,0,1) with (C,2,1,1).(C,2,2,1).

Useful Classes (for developers only)

QCMagic includes a number of useful python classes that may be used to generate or analyse data from Q-Chem calculations. Some of the more useful classes and their functions are listed below. QCManager imports most of these and thus importing QCManager sets up QCMagic for easy usage.

QCMolecule.Atom

A class used to keep track of atoms and their positions via self.Position.

QCMolecule.System

As a relatively fundamental class in QCMagic, this is used to keep track of a set of atoms contained in self.Atoms, and can be initiated with no inputs followed by addition of atoms, either individually or from a .xyz file/text. Most noteworthy are:

• numerous methods for moving and rotating both the entire system and individual atoms;
• and the method MakeInput() for generating a $molecule section of a Q-Chem input file and self.Charge, self.SpinZ to keep track of spin and charge.

QCSystemData.SystemDataSet

This class is used to keep track of multiple systems and their associated electronic states. Systems can be added either from SCF states (QCSystemData.SCFState(system,minimum)) or from a QCOutput.OutputFile object. Some important methods and variables in the class are:

• self.GeomList is an instance variable, which is a list of all systems contained in the SDS;
• self.StateList is another instance variable, which is the corresponding list of dictionaries (dicts) with each dictionary (dict) containing the minima (electronic states) associated with that system. Each state has a unique StateList index;
• RunQChem() is a class method for running a Q-Chem calculation. This takes as arguments a Queue and Job (see below), the geometry index (GI) of the System for which the calculation is to be performed, as well as a list (FindMinima) of minima to use as initial guesses given by their StateList indices. if OverwriteGeom=True, the geometry of the GI will be used, otherwise the geometry from the the next pending Job will be used.

QCOutput.OutputFile

As the name suggests, this is a class used for analysing Q-Chem output files. the Parse() method extracts data from the output file and adds it to self.Outputs. If only one calculation has been performed, self.Outputs[0] contains the data of interest. This includes steps; a list of steps corresponding to the different geometries of e.g. a geometry optimisation and the corresponding System classes can be extracted. Each step (and CurStep / Steps[-1] for a single point calculation) contains a list of Minima which can be used to instantiate a QCSystemData.SCFState object and add to an SDS. The Minima also contain orbital coefficients and distance metrics.

QCInput.Input

A class used to generate and manipulate inputs for Q-Chem calculations. An instance of this class contains a section list, self.Sections, and a section dictionary, self.SectionDict, which contain sections to be added to a Q-Chem input file. The rem section is also added as a dictionary and its elements may be changed by e.g. self.Rem['jobtype']='opt'. Sections can be added, removed or edited. self.System contains the System for which an input is to be written; n.b. this overwrites any changes to the molecule section! MakeInputText() generates the text that goes into a Q-Chem input file, but jobs may more easily be submitted directly via QCMagic.

newJob.Job

A class used to managed Q-Chem jobs. A new input can be added to the job with the AddInput() method. This input is automatically added to self.Pending from which jobs are run when the Job class instance is used to do a Q-Chem calculation.

QCQueue.QCQueue

A class used for running Q-Chem calculations. This is where input files, submission scripts, scratch directories etc. are actually written and the job submitted. This is also where details about nodes and cores are determined and may need some editing as the hardware and software used to run the calculations is changed.

QCRegisterOptions.RegisterOptions

Made to register options for scanSurface calculations (see below) but can be partly copied, imitated or imported for easy addition of own options to new scripts.

A useful way of saving and loading such class instances is to use the pickle module.