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.

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 to understand: 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 was 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.

Example command:

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

Use this to run a scan on a molecule, running qchem at points on a surface defined by the options. The above example would stretch the bond between atom0 and atom1 by 10 units in steps of 0.01, and write the results to Tag.dat.

If multiple options are used, the scan will zigzag to cover the whole surface.

Run with --help to view all possible options.

Options

--read-minimum

This is very simple to use. Just not that it is not zero-based, so --read-minimum=1 ensures only the first minimum is scanned (i.e. post-metadynamics). On the other hand --read-minimum=0 would have no effect on the scan.

--normal-rotate

This 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) your bond lengths are unusual (b) there are multiple paths between two atoms.

E.g. --normal-rotate=0,1,2,90,0.1 will:

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

--tie

An additional option --tie is soon to be included. This allows multiple parameters to change in sync with each other. For example, we might want to investigate the effect of varying two bond lengths at the same time, but don't care what happens when one is stretched and the other unchanged. This could be very useful in the context of molecular vibrations.

For each parameter you want to tie, two integers are required:

• one referring to the option;

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

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

Example command using --tie:

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

NB Ensure an equal number of steps for each tied parameter!