QCMagic

QCMagic is a set of python libraries to manipulate qchem output, orbital-coefficient etc. files. You can clone it with

  git clone ch-thom@git.csx.cam.ac.uk:QCHEM/qcmagic

The current main branch is winterClean and needs some TLC and documentation.

Documentations are underway for the side branch BCHSummer2017.

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).

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] SD_file 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 my_data.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.