pygauss 0.1.16

PYthon GAUSSian Chemical Compuation Analysis

PyGauss is designed to be an API for parsing one or more input/output files from a Gaussian quantum chemical computation and provide functionality to assess molecular geometry and electronic distribution both visually and quantitatively.

It is built on top of the cclib/chemview/chemlab suite of packages and python scientific stack and is primarily designed to be used interactively in the IPython Notebook (within which this readme was created). As shown below, a molecular optimisation can be assesed individually (much like in gaussview), but also as part of a group. The advantages of this package are then:

• Faster, more efficient analysis

• Reproducible analysis

• Trend analysis

Instillation

1. The source code is available at Github, however, the recommended way to install PyGauss is to use the Anaconda python distribution. Once downloaded a new environment can be created:

conda create -n env python=2.7

2.(L/O) If using Linux or OS X then chemlab has already been pre-built and can be installed as such:

conda install -n env -c https://conda.binstar.org/gabrielelanaro chemlab

3. PyGauss is then available for installation from PyPi after some initial dependancy installs:

conda install -n env pil
conda install -n env scipy
activate env
pip install pygauss

2.(W) Unfortuantely Windows has no pre-built installer, and so there are a few more steps to install from Github (you need to download git):

conda install -n env -c https://conda.binstar.org/gabrielelanaro cclib

conda install -n env ipython-notebook
conda install -n env numpy
conda install -n env numba
git clone https://github.com/gabrielelanaro/chemview
cd chemview
activate env
pip install .

git clone --recursive https://github.com/chemlab/chemlab.git
pip install pyopengl==3.0.2
python setup.py build_ext --inplace
add chemlab folder path to PYTHONPATH environmental variable

You should then be able to start an assessment in IPython Notebook starting with the following:

from IPython.display import display
%matplotlib inline
import pygauss as pg
folder = pg.get_test_folder()

Single Molecule Analysis

A molecule can be created containg data about the inital geometry, optimisation process and analysis of the final configuration. Molecules can be viewed statically or interactively (not currently supported by Firefox).

mol = pg.molecule.Molecule(folder,
                init_fname='CJS1_emim-cl_B_init.com',
                opt_fname=['CJS1_emim-cl_B_6-311+g-d-p-_gd3bj_opt-modredundant_difrz.log',
                           'CJS1_emim-cl_B_6-311+g-d-p-_gd3bj_opt-modredundant_difrz_err.log',
                           'CJS1_emim-cl_B_6-311+g-d-p-_gd3bj_opt-modredundant_unfrz.log'],
                freq_fname='CJS1_emim-cl_B_6-311+g-d-p-_gd3bj_freq_unfrz.log',
                nbo_fname='CJS1_emim-cl_B_6-311+g-d-p-_gd3bj_pop-nbo-full-_unfrz.log',
                alignto=[3,2,1])

#mol.show_initial(active=True)
display(mol.show_initial(zoom=0.5, rotations=[[0,0,90], [-90, 90, 0]]))
display(mol.show_optimisation(ball_stick=True, rotations=[[0,0,90], [-90, 90, 0]]))

Basic analysis of optimisation…

print('Optimised? {0}, Conformer? {1}, Energy = {2} a.u.'.format(
    mol.is_optimised(), mol.is_conformer(), round(mol.get_optimisation_E(units='hartree'),3)))
ax = mol.plot_optimisation_E(units='hartree')
ax.get_figure().set_size_inches(3, 2)

Optimised? True, Conformer? True, Energy = -805.105 a.u.

Geometric analysis…

print 'Cl optimised polar coords from aromatic ring : ({0}, {1},{2})'.format(
    *[round(i, 2) for i in mol.calc_polar_coords_from_plane(20,3,2,1)])
ax = mol.plot_opt_trajectory(20, [3,2,1])
ax.set_title('Cl optimisation path')
ax.get_figure().set_size_inches(4, 3)

Cl optimised polar coords from aromatic ring : (0.11, -116.42,-170.06)

Potential Energy Scan analysis of geometric conformers…

mol2 = pg.molecule.Molecule(folder, alignto=[3,2,1],
            pes_fname=['CJS_emim_6311_plus_d3_scan.log',
                       'CJS_emim_6311_plus_d3_scan_bck.log'])
ax = mol2.plot_pes_scans([1,4,9,10], rotation=[0,0,90], img_pos='local_maxs', zoom=0.5)
ax.set_title('Ethyl chain rotational conformer analysis')
ax.get_figure().set_size_inches(7, 3)

Natural Bond Orbital and Second Order Perturbation Theory analysis…

print '+ve charge centre polar coords from aromatic ring: ({0} {1},{2})'.format(
    *[round(i, 2) for i in mol.calc_nbo_charge_center(3, 2, 1)])
display(mol.show_nbo_charges(ball_stick=True, axis_length=0.4,
                              rotations=[[0,0,90], [-90, 90, 0]]))
display(mol.show_SOPT_bonds(min_energy=15., rotations=[[0, 0, 90]]))

+ve charge centre polar coords from aromatic ring: (0.02 -51.77,-33.15)

Multiple Computations Analysis

Multiple computations, for instance of different starting conformations, can be grouped into an Analysis class.

analysis = pg.analysis.Analysis(folder)
df, errors = analysis.add_runs(headers=['Cation', 'Anion', 'Initial'],
                               values=[['emim'], ['cl'],
                                       ['B', 'BE', 'BM', 'F', 'FE', 'FM']],
            init_pattern='CJS1_{0}-{1}_{2}_init.com',
            opt_pattern='CJS1_{0}-{1}_{2}_6-311+g-d-p-_gd3bj_opt-modredundant_unfrz.log',
            freq_pattern='CJS1_{0}-{1}_{2}_6-311+g-d-p-_gd3bj_freq_unfrz.log',
            nbo_pattern='CJS1_{0}-{1}_{2}_6-311+g-d-p-_gd3bj_pop-nbo-full-_unfrz.log')
print 'Read Errors:', errors

Read Errors: [{'Cation': 'emim', 'Initial': 'FM', 'Anion': 'cl'}]

The methods mentioned for indivdiual molecules can then be applied to all or a subset of these computations.

analysis.add_mol_property_subset('Opt', 'is_optimised', rows=[2,3])
analysis.add_mol_property('Energy (au)', 'get_optimisation_E', units='hartree')
analysis.add_mol_property('Cation chain, $\\psi$', 'calc_dihedral_angle', [1, 4, 9, 10])
analysis.add_mol_property('Cation Charge', 'calc_nbo_charge', range(1, 20))
analysis.add_mol_property('Anion Charge', 'calc_nbo_charge', [20])
analysis.add_mol_property(['Anion-Cation, $r$', 'Anion-Cation, $\\theta$', 'Anion-Cation, $\\phi$'],
                               'calc_polar_coords_from_plane', 3, 2, 1, 20)
analysis

  Anion Cation Initial   Opt  Energy (au)  Cation chain, $psi$  Cation Charge  Anion Charge  Anion-Cation, $r$  Anion-Cation, $theta$  Anion-Cation, $phi$
0    cl   emim       B   NaN     -805.105                80.794          0.888        -0.888              0.420                -123.392               172.515
1    cl   emim      BE   NaN     -805.105                80.622          0.887        -0.887              0.420                -123.449               172.806
2    cl   emim      BM  True     -805.104                73.103          0.874        -0.874              0.420                 124.121              -166.774
3    cl   emim       F  True     -805.118               147.026          0.840        -0.840              0.420                  10.393                 0.728
4    cl   emim      FE   NaN     -805.117                85.310          0.851        -0.851              0.417                 -13.254                -4.873

RadViz is a way of visualizing multi-variate data.

ax = analysis.plot_radviz_comparison('Anion', columns=range(4, 10))

The KMeans algorithm clusters data by trying to separate samples in n groups of equal variance.

kwargs = {'mtype':'optimised', 'align_to':[3,2,1],
            'rotations':[[0, 0, 90], [-90, 90, 0]],
            'axis_length':0.3}
def show_groups(df):
    for cat, gf in df.groupby('Category'):
        print 'Category {0}:'.format(cat)
        mols = analysis.yield_mol_images(rows=gf.index.tolist(), **kwargs)
        for mol, row in zip(mols, gf.index.tolist()):
            print '(row {0})'.format(row)
            display(mol)
show_groups(analysis.calc_kmean_groups('Anion', 'cl', 4, columns=range(4, 10)))

Category 0:
(row 0)

(row 1)

Category 1:
(row 3)

Category 2:
(row 2)

Category 3:
(row 4)

MORE TO COME!!