Example Assessment

After installing PyGauss you should be able to open this IPython Notebook from; https://github.com/chrisjsewell/PyGauss/blob/master/Example_Assessment.ipynb, and run the following...

from IPython.display import display, Image
%matplotlib inline
import pygauss as pg
print 'pygauss version: {}'.format(pg.__version__)

pygauss version: 0.4.3

The test folder has a number of example Gaussian outputs to play around with.

folder = pg.get_test_folder()
len(folder.list_files())

37

Note: the folder object will act identical whether using a local path or one on a server over ssh (using paramiko):

folder = pg.Folder('/path/to/folder',
                ssh_server='login.server.com',
                ssh_username='username')

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.

mol = pg.molecule.Molecule(folder_obj=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',
                atom_groups={'emim':range(20), 'cl':[20]},
                alignto=[3,2,1])

#mol.show_initial(active=True)
vdw = mol.show_initial(represent='vdw', rotations=[[0,0,90], [-90, 90, 0]])
ball_stick = mol.show_optimisation(represent='ball_stick', rotations=[[0,0,90], [-90, 90, 0]])
display(vdw, ball_stick)

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)
ax = mol.plot_freq_analysis()
ax.get_figure().set_size_inches(4, 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_obj=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(represent='ball_stick', axis_length=0.4,
                              rotations=[[0,0,90], [-90, 90, 0]]))

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

print 'H inter-bond energy = {} kJmol-1'.format(
        mol.calc_hbond_energy(eunits='kJmol-1', atom_groups=['emim', 'cl']))
print 'Other inter-bond energy = {} kJmol-1'.format(
    mol.calc_sopt_energy(eunits='kJmol-1', no_hbonds=True, atom_groups=['emim', 'cl']))
display(mol.show_sopt_bonds(min_energy=1, eunits='kJmol-1',
                            atom_groups=['emim', 'cl'],
                            no_hbonds=True,
                            rotations=[[0, 0, 90]]))
display(mol.show_hbond_analysis(cutoff_energy=5.,alpha=0.6,
                                atom_groups=['emim', 'cl'],
                                rotations=[[0, 0, 90], [90, 0, 0]]))

H inter-bond energy = 111.7128 kJmol-1
Other inter-bond energy = 11.00392 kJmol-1

Multiple Computations Analysis

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

analysis = pg.Analysis(folder_obj=folder)
errors = analysis.add_runs(headers=['Cation', 'Anion', 'Initial'],
                               values=[['emim'], ['cl'],
                                       ['B', 'BE', 'BM', 'F', 'FE']],
            init_pattern='*{0}-{1}_{2}_init.com',
            opt_pattern='*{0}-{1}_{2}_6-311+g-d-p-_gd3bj_opt*unfrz.log',
            freq_pattern='*{0}-{1}_{2}_6-311+g-d-p-_gd3bj_freq*.log',
            nbo_pattern='*{0}-{1}_{2}_6-311+g-d-p-_gd3bj_pop-nbo-full-*.log',
            alignto=[3,2,1], atom_groups={'emim':range(1,20), 'cl':[20]})

fig, caption = analysis.plot_mol_images(mtype='initial', max_cols=3,
                        info_columns=['Cation', 'Anion', 'Initial'],
                        rotations=[[0,0,90]])
print caption

Figure: (A) emim, cl, B, (B) emim, cl, BE, (C) emim, cl, BM, (D) emim, cl, F, (E) emim, cl, FE

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', 'emim')
analysis.add_mol_property('Anion Charge', 'calc_nbo_charge', 'cl')
analysis.add_mol_property(['Anion-Cation, $r$', 'Anion-Cation, $\\theta$', 'Anion-Cation, $\\phi$'],
                               'calc_polar_coords_from_plane', 3, 2, 1, 20)
analysis.add_mol_property('Anion-Cation h-bond', 'calc_hbond_energy',
                          eunits='kJmol-1', atom_groups=['emim', 'cl'])
analysis.get_table(row_index=['Anion', 'Cation', 'Initial'],
                   column_index=['Cation', 'Anion', 'Anion-Cation'])

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

NEW FEATURE: there is now an option (requiring pdflatex and ghostscript+imagemagik) to output the tables as a latex formatted image.

analysis.get_table(row_index=['Anion', 'Cation', 'Initial'],
                   column_index=['Cation', 'Anion', 'Anion-Cation'],
                   as_image=True, font_size=12)

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 into n groups of equal variance.

pg.utils.imgplot_kmean_groups(
    analysis, 'Anion', 'cl', 4, range(4, 10),
    output=['Initial'], mtype='optimised',
    rotations=[[0, 0, 90], [-90, 90, 0]],
    axis_length=0.3)

Figure: (A) BM

Figure: (A) FE

Figure: (A) B, (B) BE

Figure: (A) F

Documentation (MS Word)

After analysing the computations, it would be reasonable to want to document some of our findings. This can be achieved by outputting individual figure or table images via the folder object.

file_path = folder.save_ipyimg(vdw, 'image_of_molecule')
Image(file_path)

But you may also want to produce a more full record of your analysis, and this is where python-docx steps in. Building on this package the pygauss MSDocument class can produce a full document of your analysis.

d = pg.MSDocument()
d.add_heading('A Pygauss Example Assessment', level=1)

d.add_paragraph('We have looked at the following aspects;')
d.add_list(['geometric conformers', 'electronic structure'])

d.add_heading('Geometric Conformers', level=2)
fig, caption = analysis.plot_mol_images(max_cols=2,
                rotations=[[90,0,0], [0,0,90]],
                info_columns=['Anion', 'Cation', 'Initial'])
d.add_mpl(fig, dpi=96, height=9)
fig.clear()
d.add_markdown(caption.replace('Figure:', '**Figure:**'))
d.add_paragraph()
df = analysis.get_table(columns=['Anion Charge', 'Cation Charge',
                                 'Energy (au)'],
                   row_index=['Anion', 'Cation', 'Initial'])
d.add_dataframe(df, incl_indx=True, style='Medium Shading 1 Accent 1')
d.add_markdown('**Table:** Analysis of Conformer Charge')

d.save('exmpl_assess.docx')

Which gives us the following:

MORE TO COME!!