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

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!