pygauss 0.1.14

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/chemlab/chemview suite of packages and python scientific stack and is primarily designed to be used interactively in the IPython Notebook (within which this readme has been written). 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

Detail instillation…

pip install pygauss

conda install -c http://conda.binstar.org/gabrielelanaro chemlab

You should then be able to start an ipython notebook…

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.

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

a

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'}]

analysis.add_mol_property('Opt', 'is_optimised')
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  True     -805.105                80.794          0.888        -0.888              0.420                -123.392               172.515
1    cl   emim      BE  True     -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  True     -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 3)

Category 1:
(row 0)

(row 1)

Category 2:
(row 2)

Category 3:
(row 4)