Skip to main content

PYthon GAUSSian DFT output analysis

Project description

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]]))
https://github.com/chrisjsewell/PyGauss/readmeraw/master/output_6_0.png https://github.com/chrisjsewell/PyGauss/readmeraw/master/output_6_1.png

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.
https://github.com/chrisjsewell/PyGauss/readmeraw/master/output_8_1.png

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)
https://github.com/chrisjsewell/PyGauss/readmeraw/master/output_10_1.png

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)
https://github.com/chrisjsewell/PyGauss/readmeraw/master/output_12_0.png

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)
https://github.com/chrisjsewell/PyGauss/readmeraw/master/output_14_1.png https://github.com/chrisjsewell/PyGauss/readmeraw/master/output_14_2.png

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))
https://github.com/chrisjsewell/PyGauss/readmeraw/master/output_20_0.png

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)
https://github.com/chrisjsewell/PyGauss/readmeraw/master/output_22_1.png
Category 1:
(row 0)
https://github.com/chrisjsewell/PyGauss/readmeraw/master/output_22_3.png
(row 1)
https://github.com/chrisjsewell/PyGauss/readmeraw/master/output_22_5.png
Category 2:
(row 2)
https://github.com/chrisjsewell/PyGauss/readmeraw/master/output_22_7.png
Category 3:
(row 4)
https://github.com/chrisjsewell/PyGauss/readmeraw/master/output_22_9.png

Project details


Download files

Download the file for your platform. If you're not sure which to choose, learn more about installing packages.

Source Distribution

pygauss-0.1.8.zip (5.0 MB view hashes)

Uploaded Source

Supported by

AWS AWS Cloud computing and Security Sponsor Datadog Datadog Monitoring Fastly Fastly CDN Google Google Download Analytics Microsoft Microsoft PSF Sponsor Pingdom Pingdom Monitoring Sentry Sentry Error logging StatusPage StatusPage Status page