In [252]:
import pygauss as pg
import pandas as pd
import matplotlib.pyplot as plt
from IPython.display import display
%matplotlib inline
import warnings
warnings.filterwarnings("ignore")
print 'pygauss version: {}'.format(pg.__version__)
In [305]:
import os, platform
home = os.path.expanduser("~")
if platform.system() == 'Windows':
pg.set_imagik_exe('convert_pdf')
onedrive = 'SkyDrive'
else:
onedrive = 'OneDrive'
inpath = os.path.join(home, onedrive,
'Imperial_2014-15', 'Project', 'Gaussian_files')
In [444]:
il_coronene = pg.Analysis(os.path.join(inpath, 'Coronene-IL'))
cerrs = il_coronene.add_runs(
headers=['Surface', 'Cation', 'Anion', 'Diffuse', 'Initial'],
values=[['none'], ['emim'], ['none'], ['g', '+g'], ['none']],
init_pattern='*_{1}_-_init.com',
opt_pattern='*_{1}_-_6-311{3}-*_opt_.log',
freq_pattern='*_{1}_-_6-311{3}-*_freq_.log',
#nbo_pattern='*_{1}_-_6-311{3}-*_pop-nbo-full-_.log',
atom_groups={'emim':range(1, 20)}, alignto=[3,2,1])
aerrs = il_coronene.add_runs(
headers=['Surface', 'Cation', 'Anion', 'Diffuse', 'Initial'],
values=[['none'], ['none'], ['cl', 'bf', 'etso'], ['g', '+g'], ['none']],
init_pattern='*_{2}_-_init.com',
opt_pattern='*_{2}_-_6-311{3}-d-p-_gd3bj_opt_.log',
freq_pattern='*_{2}_-_6-311{3}-d-p-_gd3bj_freq_.log')
serrs = il_coronene.add_runs(
headers=['Surface', 'Cation', 'Anion', 'Diffuse', 'Initial'],
values=[['coro'], ['none'], ['none'], ['g', '+g'], ['none']],
init_pattern='*_{0}_-_init.com',
opt_pattern='*_{0}_-_6-311{3}-d-p-_gd3bj_opt_*.log',
freq_pattern='*_{0}_-_6-311{3}-d-p-_gd3bj_freq_*.log')
error_df = pd.concat([cerrs, aerrs, serrs])
for anion, emim_atoms in zip(['cl', 'bf'],
[range(38, 57), range(42, 61)]):
atom_groups={'coro':range(1,37), 'coro_plane':[1,12,28], 'emim':emim_atoms, 'ce':range(1,37)+emim_atoms,
'emim_ring':list(reversed(emim_atoms[0:3])), 'c2':[emim_atoms[2]],
'cl':[37], 'bf':range(37, 42), 'coro_plane':[4,7,9]}
emim_hs = [(2,7), (4,6), (5,5), ('6a',17), ('6b',18), ('6c',19),
('7a',11), ('7b',12), ('8a',13), ('8b',14), ('8c',15)]
for name, number in emim_hs:
atom_groups['h{}'.format(name)] = emim_atoms[number-1]
label = 'H{}_charge'.format(name)
errors = il_coronene.add_runs(
headers=['Surface', 'Cation', 'Anion', 'Initial', 'Diffuse'],
values=[['coro'], ['emim'], [anion],
['P-'+str(a) for a in range(0,360,30)]+['O-O', 'O-S1'],
['+g']],
init_pattern='CJS2_{0}-{1}-{2}_{3}_init.com',
opt_pattern='CJS2_{0}-{1}-{2}_{3}_6-311{4}-d-p-_gd3bj_opt_*.log',
freq_pattern='CJS2_{0}-{1}-{2}_{3}_6-311{4}-d-p-_gd3bj_freq_*.log',
nbo_pattern='CJS6_{0}-{1}-{2}_{3}_6-311{4}-d-p-_gd3bj_pop-nbo-full-_*.log',
atom_groups=atom_groups,
alignto=[4,7,9], add_if_error=True, ipython_print=True)
error_df = pd.concat([error_df, errors])
atom_groups={'coro':range(1,37), 'coro_plane':[1,12,28], 'emim':range(37, 56),
'ce':range(1,37)+range(37, 56), 'emim_ring':[39,38,37],
'cl':[56], 'bf':range(56, 61), 'coro_plane':[4,7,9], 'c2':[39]}
emim_hs = [(2,7), (4,6), (5,5), ('6a',17), ('6b',18), ('6c',19),
('7a',11), ('7b',12), ('8a',13), ('8b',14), ('8c',15)]
for name, number in emim_hs:
atom_groups['h{}'.format(name)] = 36+number
errors = il_coronene.add_runs(
headers=['Surface', 'Cation', 'Anion', 'Initial', 'Diffuse'],
values=[['coro'], ['emim'], ['cl', 'bf'],
['O-S2'], ['+g']],
init_pattern='CJS2_{0}-{1}-{2}_{3}_init.com',
opt_pattern='CJS2_{0}-{1}-{2}_{3}_6-311{4}-d-p-_gd3bj_opt_*.log',
freq_pattern='CJS2_{0}-{1}-{2}_{3}_6-311{4}-d-p-_gd3bj_freq_*.log',
nbo_pattern='CJS6_{0}-{1}-{2}_{3}_6-311{4}-d-p-_gd3bj_pop-nbo-full-_*.log',
atom_groups=atom_groups,
alignto=[4,7,9], add_if_error=True, ipython_print=True)
error_df = pd.concat([error_df, errors])
for anion, emim_atoms in zip(['etso'],
[range(49,68)]):
atom_groups={'coro':range(1,37), 'coro_plane':[1,12,28], 'emim':emim_atoms, 'ce':range(1,37)+emim_atoms, 'emim_ring':list(reversed(emim_atoms[0:3])),
'etso':range(37,49), 'coro_plane':[4,7,9], 'c2':[emim_atoms[2]]}
emim_hs = [(2,7), (4,6), (5,5), ('6a',17), ('6b',18), ('6c',19),
('7a',11), ('7b',12), ('8a',13), ('8b',14), ('8c',15)]
for name, number in emim_hs:
atom_groups['h{}'.format(name)] = emim_atoms[number-1]
errors = il_coronene.add_runs(
headers=['Surface', 'Cation', 'Anion', 'Initial', 'Diffuse'],
values=[['coro'], ['emim'], [anion],
['P2-'+str(a) for a in range(0,360,30)],
['+g']],
init_pattern='CJS2_{0}-{1}-{2}_{3}_init.com',
opt_pattern='CJS2_{0}-{1}-{2}_{3}_6-311{4}-d-p-_gd3bj_opt_*.log',
freq_pattern='CJS2_{0}-{1}-{2}_{3}_6-311{4}-d-p-_gd3bj_freq_*.log',
nbo_pattern='CJS6_{0}-{1}-{2}_{3}_6-311{4}-d-p-_*_pop-nbo-full-_*.log',
atom_groups=atom_groups,
alignto=[4,7,9], add_if_error=True, ipython_print=True)
error_df = pd.concat([error_df, errors])
print 'Read Errors:'
error_df[error_df.Type!='nbo']
Out[444]:
In [445]:
il_coronene._df['InPlane'] = il_coronene._df.Initial.str.contains('P')
il_coronene.get_table(filters={'Anion':'cl', 'Surface':'coro'})
Out[445]:
In [446]:
fig, caption = il_coronene.plot_mol_images(mtype='highlight-initial',
highlight=['emim', 'cl'], alpha=0.1, transparent=False, info_incl_id=True,
info_columns=['Initial'], max_cols=3, padding=(1,1),
filters={'Surface':'coro', 'Anion':'cl', 'Diffuse':'+g', 'InPlane':False},
rotations=[[0, 90, 90]], label_size=0, align_to='emim_ring', width=1000)
print caption
fig.savefig('out_of_plane_anion_positions.png', dpi=200)
In [447]:
il_coro_confs = il_coronene.copy()
non_conformers = il_coro_confs.remove_non_conformers()
#add cl back because theres a problem with it being a single atom
#and regestering it as optimised/passing frequency analysis
il_coro_confs.add_runs(
headers=['Surface', 'Anion', 'Cation', 'Diffuse', 'Initial'],
values=[['none'], ['cl'], ['none'], ['g', '+g'], ['none']],
init_pattern='*_{1}_-_init.com',
opt_pattern='*_{1}_-_6-311{3}-d-p-_gd3bj_opt_.log',
freq_pattern='*_{1}_-_6-311{3}-d-p-_gd3bj_freq_.log')
#similarly with coronene and a diffuse function
il_coro_confs.add_runs(
headers=['Surface', 'Anion', 'Cation', 'Diffuse', 'Initial'],
values=[['coro'], ['none'], ['none'], ['+g'], ['none']],
init_pattern='*_{0}_-_init.com',
opt_pattern='*_{0}_-_6-311{3}-d-p-_gd3bj_opt_*.log',
freq_pattern='*_{0}_-_6-311{3}-d-p-_gd3bj_freq_*.log')
Out[447]:
In [448]:
for anion in ['cl', 'bf', 'etso']:
s_id, c_id, a_id = il_coro_confs.get_ids(['Surface', 'Anion', 'Cation', 'Diffuse'],
[['coro', 'none', 'none', '+g'],
['none', anion, 'none', '+g'], ['none', 'none', 'emim', '+g']])
il_coro_confs.add_mol_property_subset(
'Energy of Association (kJmol^{-1}) Uncorrected',
'get_opt_energy',
kwargs={'units':'kJmol-1'},
filters={'Surface':'coro', 'Diffuse':'+g',
'Cation':'emim', 'Anion':anion},
relative_to_rows=[s_id, c_id, a_id])
il_coro_confs.add_mol_property_subset(
'Energy of Association (kJmol^{-1}) Corrected',
'get_opt_energy',
kwargs={'units':'kJmol-1', 'zpe_correct':True},
filters={'Surface':'coro', 'Diffuse':'+g',
'Cation':'emim', 'Anion':anion},
relative_to_rows=[s_id, c_id, a_id])
for diff in ['+g', 'g']:
ids = il_coro_confs.get_ids(
['Surface', 'Anion', 'Cation', 'Diffuse', 'Initial'],
[['coro', 'none', 'none', diff, 'none'],
['none', 'none', 'emim', diff, 'none'],
['none', 'cl', 'none', diff, 'none'],
['none', 'bf', 'none', diff, 'none'],
['none', 'etso', 'none', diff, 'none']])
df = il_coro_confs.remove_rows(ids)
In [449]:
df = il_coro_confs.get_table(
precision=3,
column_index=['Energy of Association (kJmol^{-1})'])
df.sort([(' ', 'Anion'), ('Energy of Association (kJmol^{-1})','Corrected')])
Out[449]:
In [450]:
for anion in ['cl']:
fig, caption = il_coro_confs.plot_mol_images(mtype='highlight',
highlight=['emim', anion], alpha=0.1, transparent=True, info_incl_id=True,
info_columns=['Initial'], max_cols=3, padding=(1,1),
filters={'Surface':'coro', 'Anion':anion, 'Diffuse':'+g'},
rotations=[[0, 0, 90], [0, 90, 90]], label_size=8, align_to='emim_ring')
print caption
In [451]:
fig, caption = il_coro_confs.plot_mol_images(mtype='highlight',
highlight=['emim', anion], alpha=0.1, transparent=True, info_incl_id=True,
info_columns=['Initial'], max_cols=3, padding=(1,1),
rows=[10,11,12],
rotations=[[0, 0, 90], [0, 90, 90]], label_size=8, align_to='emim_ring')
print caption
In [452]:
fig, caption = il_coro_confs.plot_mol_images(mtype='highlight',
highlight=['emim', anion], alpha=0.1, transparent=True, info_incl_id=True,
info_columns=['Initial'], max_cols=3, padding=(1,1),
rows=[16,17],
rotations=[[0, 0, 90], [0, 90, 90]], label_size=8, align_to='emim_ring')
print caption
In [453]:
fig, caption = il_coro_confs.plot_mol_images(mtype='highlight',
highlight=['emim', anion], alpha=0.1, transparent=True, info_incl_id=True,
info_columns=['Initial'], max_cols=3, padding=(1,1),
rows=[23],
rotations=[[0, 0, 90], [0, 90, 90]], label_size=8, align_to='emim_ring')
print caption
In [454]:
fig, caption = il_coro_confs.plot_mol_images(mtype='highlight',
highlight=['emim', anion], alpha=0.1, transparent=True, info_incl_id=True,
info_columns=['Initial'], max_cols=3, padding=(1,1),
rows=[13,14,15],
rotations=[[0, 0, 90], [0, 90, 90]], label_size=8, align_to='emim_ring')
print caption
In [455]:
fig, caption = il_coro_confs.plot_mol_images(mtype='highlight',
highlight=['emim', anion], alpha=0.1, transparent=True, info_incl_id=True,
info_columns=['Initial'], max_cols=3, padding=(1,1),
rows=[19,20,21],
rotations=[[0, 0, 90], [0, 90, 90]], label_size=8, align_to='emim_ring')
print caption
In [456]:
fig, caption = il_coro_confs.plot_mol_images(mtype='highlight',
highlight=['emim', anion], alpha=0.1, transparent=True, info_incl_id=True,
info_columns=['Initial'], max_cols=3, padding=(1,1),
rows=[22],
rotations=[[0, 0, 90], [0, 90, 90]], label_size=8, align_to='emim_ring')
print caption
final conformers
In [457]:
fig, caption = il_coro_confs.plot_mol_images(mtype='highlight',
highlight=['emim', anion], alpha=0.1, transparent=True, info_incl_id=True,
info_columns=['Initial'], max_cols=3, padding=(1,1),
rows=[10,16,23,13,19,22],
rotations=[[0, 0, 90], [0, 90, 90]], label_size=8, align_to='emim_ring')
print caption
In [458]:
for anion in ['bf']:
fig, caption = il_coro_confs.plot_mol_images(mtype='highlight',
highlight=['emim', anion], alpha=0.1, transparent=True, info_incl_id=True,
info_columns=['Initial'], max_cols=3, padding=(1,1),
filters={'Surface':'coro', 'Anion':anion, 'Diffuse':'+g'},
rotations=[[0, 0, 90], [0, 90, 90]], label_size=8, align_to='emim_ring')
print caption
final conformers
In [459]:
fig, caption = il_coro_confs.plot_mol_images(mtype='highlight',
highlight=['emim', 'bf'], alpha=0.1, transparent=True, info_incl_id=True,
info_columns=['Initial'], max_cols=3, padding=(1,1),
rows=[24,28,31,35],
rotations=[[0, 0, 90], [0, 90, 90]], label_size=8, align_to='emim_ring')
print caption
In [460]:
for anion in ['etso']:
fig, caption = il_coro_confs.plot_mol_images(mtype='highlight',
highlight=['emim', anion], alpha=0.1, transparent=True, info_incl_id=True,
info_columns=['Initial'], max_cols=2, padding=(1,1),
filters={'Surface':'coro', 'Anion':anion, 'Diffuse':'+g'},
rotations=[[0, 0, 90], [0, 90, 90]], label_size=8, align_to='emim_ring')
print caption
In [461]:
fig, caption = il_coro_confs.plot_mol_images(mtype='highlight',
highlight=['emim', 'etso'], alpha=0.1, transparent=True, info_incl_id=True,
info_columns=['Initial'], max_cols=3, padding=(1,1),
rows=[43, 44, 46, 48, 49, 50],
rotations=[[0, 0, 90], [0, 90, 90]], label_size=8, align_to='emim_ring')
print caption
NB: 43 is marginally different to 44 in etso4 chain orientation
In [462]:
#leaving out 22 where cl below plane
df = il_coro_confs.remove_rows(
set(il_coro_confs.get_table().index).difference(
[10,16,23,13,19,24,28,31,35,43,46,48,49,50]))
In [463]:
def il_coro_show(anion, letter, mtype='highlight', padding=(1,2), label_size=10):
fig, caption = il_coro_confs.plot_mol_images(mtype=mtype, letter_prefix=letter,
highlight=['emim', anion], alpha=0.2, transparent=False, info_incl_id=True,
info_columns=['Initial'], max_cols=2, padding=padding,
filters={'Surface':'coro', 'Anion':anion, 'Diffuse':'+g'},
rotations=[[0, 0, 90], [0, 90, 90]], label_size=label_size, align_to='emim_ring',
sort_columns=['Energy of Association (kJmol^{-1}) Corrected'], width=1000)
print caption
return fig
fig_ecl = il_coro_show('cl', '9')
fig_ebf = il_coro_show('bf','10')
fig_eetso = il_coro_show('etso', '11', padding=(-10,2), label_size=8)
In [464]:
def il_coro_show2(anion, letter, mtype='hbond', padding=(1,2), label_size=10):
fig, caption = il_coro_confs.plot_mol_images(mtype=mtype, letter_prefix=letter,sopt_min_energy=1,
atom_groups=['ce', anion],
highlight=['emim', anion], alpha=0.6, transparent=False, info_incl_id=True,
info_columns=['Initial'], max_cols=3, padding=padding,
filters={'Surface':'coro', 'Anion':anion, 'Diffuse':'+g'},
rotations=[[80, 0, 0]], label_size=label_size, align_to='emim_ring',
sort_columns=['Energy of Association (kJmol^{-1}) Corrected'], width=1000)
print caption
return fig
fig_ecl = il_coro_show2('cl', '9')
fig_ebf = il_coro_show2('bf','10')
fig_eetso = il_coro_show2('etso', '11', padding=(1,2), label_size=10)
In [465]:
fig_ccl = il_coro_show('cl', '9', mtype='nbo')
fig_cbf = il_coro_show('bf','10', mtype='nbo')
fig_cetso = il_coro_show('etso', '11', mtype='nbo')
In [466]:
import string
il_coro_confs._df['Order'] = il_coro_confs._df['Anion'].replace(['cl', 'bf', 'etso'], [1,2,3])
il_coro_confs._df.sort(['Order', 'Energy of Association (kJmol^{-1}) Corrected'], inplace=True)
il_coro_confs._df['Letter']=['9A', '9B', '9C', '9D', '9E',
'10A', '10B', '10C', '10D', '11A', '11B', '11C', '11D', '11E']
il_coro_confs._df['Initials']=['P270,P300,P330', 'O-S1', 'P180,P210', 'P90,P120,P150', 'P0,P30,P60',
'P330', 'P210,P240', 'P120', 'P0,P60',
'P270', 'P240', 'P300', 'P180', 'P90,P120,P150']
il_coro_confs.get_table()
Out[466]:
In [467]:
import string
import copy
def plt_value(analysis, x_col, y_col,
ytick_names, ytick_labels,
x_label=None, y_label=None,
filters={}, ax=None, combine_gap=1.):
df = analysis.get_table(filters=filters).sort([y_col, x_col])
df = df.replace(ytick_names, range(1,len(ytick_names)+1))
df = df.sort(y_col)
ax = df.plot(kind='scatter', y=y_col, x=x_col, ax=ax)
ax.grid(True)
ax.set_yticks(range(1,len(ytick_names)+1))
ax.set_yticklabels(ytick_labels)
ax.set_ylabel('')
if x_label: ax.set_xlabel(x_label)
X0, Y0 = None, None
for X, Y, Z in zip(df[x_col], df[y_col], df.Letter):
# Annotate the points 5 _points_ above and to the left of the vertex
if not X0==None:
if Y==Y0 and X-X0<combine_gap:
A0.set_text(','.join([A0.get_text(),Z]))
X0, Y0 = X,Y
continue
A0 = ax.annotate(Z, xy=(X,Y), xytext=(-5, 5), ha='right',
textcoords='offset points')
X0, Y0 = X,Y
ax.invert_yaxis()
return ax
ax_energy = plt_value(il_coro_confs, 'Energy of Association (kJmol^{-1}) Corrected', 'Anion',
['cl', 'bf', 'etso'],
['$[C_{24}H_{12}][C2C1im]^+[Cl]^-$', '$[C_{24}H_{12}][C2C1im]^+[BF_4]^-$',
'$[C_{24}H_{12}][C2C1im]^+[EtSO_4]^-$'],
x_label='Energy of Association ($kJmol^{-1}$)',
filters={'Cation':'emim'}, combine_gap=3)
In [468]:
fig1,(ax1,ax2) = plt.subplots(1, 2, sharey=True, gridspec_kw={'width_ratios':[7, 2]})
ax1 = plt_value(il_coro_confs, 'Energy of Association (kJmol^{-1}) Corrected', 'Anion',
['cl', 'bf', 'etso'], ['$Cl$', '$BF_4$', '$EtSO_4$'],
x_label='Energy of Association ($kJmol^{-1}$)',
filters={'Cation':'emim'}, ax=ax1)
ax2 = plt_value(il_coro_confs, 'Energy of Association (kJmol^{-1}) Corrected', 'Anion',
['cl', 'bf', 'etso'], ['$Cl$', '$BF_4$', '$EtSO_4$'],
x_label='Energy of Association ($kJmol^{-1}$)',
filters={'Cation':'emim'}, ax=ax2)
ax1.set_ylim(0.5, 3.5)
ax1.set_xlim(-490, -420)
ax2.set_xlim(-350, -330)
ax2.set_xticks([-350, -340, -330])
ax2.set_ylabel('')
ax2.set_xlabel('')
ax1.invert_yaxis()
# hide the spines between ax and ax2
ax1.spines['right'].set_visible(False)
ax2.spines['left'].set_visible(False)
ax1.yaxis.tick_left()
ax1.tick_params(labeltop='off') # don't put tick labels at the top
#ax2.yaxis.tick_right()
# Make the spacing between the two axes a bit smaller
plt.subplots_adjust(wspace=0.2)
d = .015 # how big to make the diagonal lines in axes coordinates
# arguments to pass plot, just so we don't keep repeating them
kwargs = dict(transform=ax1.transAxes, color='k', clip_on=False)
ax1.plot((1-d,1+d),(-d,+d), **kwargs) # top-left diagonal
ax1.plot((1-d,1+d),(1-d,1+d), **kwargs) # bottom-left diagonal
kwargs.update(transform=ax2.transAxes) # switch to the bottom axes
ax2.plot((-d*7/2.,d*7/2.),(-d,+d), **kwargs) # top-right diagonal
ax2.plot((-d*7/2.,d*7/2.),(1-d,1+d), **kwargs) # bottom-right diagonal
Out[468]:
In [469]:
il_coro_confs.add_mol_property_subset('cat_coro_dist', 'calc_min_dist',
args=['emim', 'coro'], kwargs={'units':'Angstrom'})
il_coro_confs.add_mol_property_subset('ring_coro_dist', 'calc_min_dist',
args=['emim_ring', 'coro'], kwargs={'units':'Angstrom'})
for ana in ['cl', 'bf', 'etso']:
il_coro_confs.add_mol_property_subset('an_coro_dist', 'calc_min_dist',
args=[ana, 'coro'], kwargs={'units':'Angstrom'},
filters={'Anion':ana})
il_coro_confs.add_mol_property_subset('an_cat_dist', 'calc_min_dist',
args=[ana, 'emim'], kwargs={'units':'Angstrom'},
filters={'Anion':ana})
il_coro_confs.get_table()
Out[469]:
In [535]:
fig_dist,(ax1,ax2,ax3) = plt.subplots(1, 3, sharey=True, sharex=True)
df = il_coro_confs.get_table()
ax1.scatter(df['cat_coro_dist'], range(1,15))
ax1.set_xlabel('[C2C1im]$^+$ to [C$_{24}$H$_{12}$]')
ax1.grid(True)
ax2.scatter(df['an_coro_dist'], range(1,15))
ax2.set_xlabel('Anion to [C$_{24}$H$_{12}$]')
ax2.grid(True)
ax3.scatter(df['an_cat_dist'], range(1,15))
ax3.set_xlabel('[C2C1im]$^+$ to Anion')
ax3.grid(True)
ax1.set_yticks(range(1,15))
ax1.set_xticks([1.6, 2, 2.4, 2.8, 3.2])
ion_dict={
'emim':'[C2C1im]$^+$',
'omim':'[C4C1im]$^+$',
'cl':'[Cl]$^-$',
'bf':'[BF$_4$]$^-$',
'etso':'[C$_2$SO$_4$]$^-$'}
ylabels=[]
last=None
for ix, rw in df.iterrows():
if (rw.Cation, rw.Anion)!=last:
ylabels.append(''.join([ion_dict[rw.Anion],' ',rw.Letter]))
# ylabels.append(''.join([ion_dict[rw.Cation],ion_dict[rw.Anion],' ',rw.Letter]))
else:
ylabels.append(rw.Letter)
last = (rw.Cation, rw.Anion)
ax1.set_yticklabels(ylabels)
# The big subplot for x label
ax = fig_dist.add_subplot(111)
ax.tick_params(top='off', bottom='off', left='off', right='off',
labelbottom='on', labelleft='on', pad=35)
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.set_frame_on(False)
ax.set_xlabel(r'Minimum Distance ($\AA$)')
ax1.invert_yaxis()
In [553]:
fig_dist,ax1 = plt.subplots(1, 1, sharey=True, sharex=True)
df = il_coro_confs.get_table()
ax1.scatter(df['cat_coro_dist'], range(1,15))
ax1.plot(df['cat_coro_dist'], range(1,15), '--')
ax1.scatter(df['an_coro_dist'], range(1,15), c='r', marker='^')
ax1.plot(df['an_coro_dist'], range(1,15), 'r--')
ax1.scatter(df['an_cat_dist'], range(1,15), c='g', marker='s')
ax1.plot(df['an_cat_dist'], range(1,15), 'g--')
ax1.set_xlabel(r'Minimum Distance ($\AA$)')
ax1.grid(True)
ax1.set_yticks(range(1,15))
ax1.set_ybound(0.5,15)
ax1.set_xticks([1.6, 2, 2.4, 2.8, 3.2])
ion_dict={
'emim':'[C2C1im]$^+$',
'omim':'[C4C1im]$^+$',
'cl':'[Cl]$^-$',
'bf':'[BF$_4$]$^-$',
'etso':'[C$_2$SO$_4$]$^-$'}
ylabels=[]
last=None
for ix, rw in df.iterrows():
if (rw.Cation, rw.Anion)!=last:
ylabels.append(''.join([ion_dict[rw.Anion],' ',rw.Letter]))
# ylabels.append(''.join([ion_dict[rw.Cation],ion_dict[rw.Anion],' ',rw.Letter]))
else:
ylabels.append(rw.Letter)
last = (rw.Cation, rw.Anion)
ax1.set_yticklabels(ylabels)
ax1.legend(['Cation-Coronene', 'Anion-Coronene', 'Cation-Anion'], bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.,
prop={'size':12})
ax1.invert_yaxis()
fig_dist.set_size_inches(3,4)
In [488]:
il_coro_confs.add_mol_property_subset('coro_charge', 'calc_nbo_charge',
args=['coro'])
il_coro_confs.add_mol_property_subset('cat_charge', 'calc_nbo_charge',
args=['emim'])
il_coro_confs.add_mol_property_subset('C2_charge', 'calc_nbo_charge',
args=['c2'])
emim_hs = [(2,7), (4,6), (5,5), ('6a',17), ('6b',18), ('6c',19),
('7a',11), ('7b',12), ('8a',13), ('8b',14), ('8c',15)]
for name, number in emim_hs:
il_coro_confs.add_mol_property_subset('H{}_charge'.format(name), 'calc_nbo_charge',
args=['h{}'.format(name)])
for ana in ['cl', 'bf', 'etso']:
il_coro_confs.add_mol_property_subset('an_charge', 'calc_nbo_charge',
args=[ana], filters={'Anion':ana})
il_coro_confs._df['Charge_Transfer'] = 1-il_coro_confs._df['cat_charge']
il_coro_confs.get_table(precision=3)
Out[488]:
In [472]:
ax_charge = plt_value(il_coro_confs, 'Charge_Transfer', 'Anion',
['cl', 'bf', 'etso'],
['$[C_{24}H_{12}][C2C1im]^+[Cl]^-$', '$[C_{24}H_{12}][C2C1im]^+[BF_4]^-$',
'$[C_{24}H_{12}][C2C1im]^+[EtSO_4]^-$'],
x_label=r'Charge Transfer to Cation ($\Delta q_{cation}$)',
filters={'Cation':'emim'}, combine_gap=0.002)
In [473]:
il_coro_confs.add_mol_property_subset('H-bond coro-cat', 'calc_hbond_energy', args=[['coro', 'emim']])
il_coro_confs.add_mol_property_subset('Other-bond coro-cat', 'calc_sopt_energy', args=[['coro', 'emim']],
kwargs={'no_hbonds':True})
Out[473]:
In [474]:
for ani in ['cl', 'bf', 'etso']:
il_coro_confs.add_mol_property_subset('H-bond coro-an', 'calc_hbond_energy', args=[['coro', ani]],
filters={'Anion':ani})
il_coro_confs.add_mol_property_subset('Other-bond coro-an', 'calc_sopt_energy', args=[['coro', ani]],
kwargs={'no_hbonds':True}, filters={'Anion':ani})
il_coro_confs.add_mol_property_subset('H-bond cat-an', 'calc_hbond_energy', args=[['emim', ani]],
filters={'Anion':ani})
il_coro_confs.add_mol_property_subset('Other-bond cat-an', 'calc_sopt_energy', args=[['emim', ani]],
kwargs={'no_hbonds':True}, filters={'Anion':ani})
il_coro_confs.get_table(precision=3)
Out[474]:
In [475]:
for anion, mol in zip(il_coro_confs._df.Anion, il_coro_confs._df.Molecule):
print '-----'
a= mol.get_hbond_analysis(atom_groups=['coro', anion])
a['D']=[str(_) for _ in a.Donors]
a['A']=[str(_) for _ in a.Acceptors]
print a.groupby(['D', 'A']).E2.sum()
In [525]:
fig_sopt,(ax1,ax2) = plt.subplots(1, 2, sharey=True, gridspec_kw={'width_ratios':[1, 2]})
df = il_coro_confs._df.sort(['Order', 'Letter'])
df.plot(kind='barh',
y=['Other-bond {}'.format(t) for t in ['cat-an', 'coro-cat', 'coro-an']] , x='Letter',
ax=ax1, legend=False, fontsize=12)
for container in ax1.containers:
plt.setp(container, height=0.2)
ax1.grid(True)
ax1.set_xlim(0, 70)
ax1.set_xticks([0,20,40,60])
ax1.set_ylabel('')
ax1.set_xlabel('Other', fontdict={'fontsize':11})
ax1.invert_xaxis()
df.plot(kind='barh', y=['H-bond {}'.format(t) for t in ['cat-an', 'coro-cat', 'coro-an']], x='Letter',
ax=ax2, legend=False, fontsize=12)
for container in ax2.containers:
plt.setp(container, height=0.2)
ax2.grid(True)
ax2.set_xlim(0, 140)
ax2.set_ylabel('')
ax2.set_xlabel('H-Bond', fontdict={'fontsize':11})
ax2.legend(['Cation-Anion', 'Cation-Coronene', 'Anion-Coronene'], bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.,
prop={'size':12})
ion_dict={
'emim':'[C2C1im]^+',
'omim':'[C4C1im]^+',
'cl':'[Cl]^-',
'bf':'[BF_4]^-',
'etso':'[EtSO_4]^-'}
ylabels=[]
last=None
for ix, rw in df.iterrows():
if (rw.Cation, rw.Anion)!=last:
ylabels.append(''.join(['$',ion_dict[rw.Cation],ion_dict[rw.Anion],'$ ',rw.Letter]))
else:
ylabels.append(rw.Letter)
last = (rw.Cation, rw.Anion)
ax1.set_yticklabels(ylabels)
ax1.invert_yaxis()
# The big subplot for x label
ax = fig_sopt.add_subplot(111)
ax.tick_params(top='off', bottom='off', left='off', right='off',
labelbottom='on', labelleft='on', pad=30)
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.set_frame_on(False)
ax.set_xlabel('Inter-Molecular Donor-Acceptor Interactions ($kJmol^{-1}$)', fontdict={'fontsize':12}) #from Second-Order Perturbative
fig_sopt.set_size_inches(6,4)
fig_sopt.savefig('test', dpi=400, bbox_inches='tight')
In [477]:
for anion, l in zip(['cl', 'bf', 'etso'], ['A', 'F', 'J'], []):
fig, caption = il_coro_confs.plot_mol_graphs(gtype='dos', max_cols=3,
per_energy=1, lbound=-25, ubound=5, start_letter=l,
filters={'Anion':anion}, legend_size=8,
info_incl_id=True,
atom_groups=['emim', anion], group_colors=['blue', 'magenta'],
group_labels=['EMIM', anion.upper()], group_fill=False)
for ax in fig.get_axes():
l = ax.get_legend()
if l:
l.set_visible(False)
else:
ax0.legend(['HOMO', 'LUMO', '[C2C1im]', anion.upper()], bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.,
prop={'size':10})
break
ax0 = ax
In [478]:
il_coro_confs.get_table(filters={'Cation':'emim'}).sort(['Order', 'Letter'])
Out[478]:
In [479]:
data = il_coro_confs.get_table(filters={'Cation':'emim'}).sort(['Order', 'Letter'])
ilg_df=pd.DataFrame()
ilg_df['Conformer'] = data.Letter
ilg_df['Starting Positions of Anion'] = data.Initials
energies, hbond, dist = [], [], []
h_dict_ca_hbond = {'cl':18.159,'bf':20.292,'etso':36.861}
h_dict_ca_dist = {'cl':2.697,'bf':2.211,'etso':2.188}
for anion, df in data.groupby('Anion', sort=False):
e = df['Energy of Association (kJmol^{-1}) Corrected']
e = e - e.min()
energies.extend(e.tolist())
h=df['H-bond coro-an'] - h_dict_ca_hbond[anion]
hbond.extend(h.tolist())
d=df['an_coro_dist'] - h_dict_ca_dist[anion]
dist.extend(d.tolist())
ilg_df[r'Energy $\Delta E$ (kJmol^{-1})']=energies
ilg_df[r'Minimum Distance ($Angstrom$) ring-cat delta'] = data.ring_coro_dist - 3.40985687794
ilg_df[r'Minimum Distance ($Angstrom$) coro_an delta'] = dist
ilg_df['coro-an H delta']=hbond
ilg_df[r'C2'] = data.C2_charge
ilg_df[r'H2'] = data.H2_charge
ilg_df['coro charge'] = data.coro_charge
#ilg_df[r'Charge Transfer to coronene ($\Delta$ q_{cation})'] = data.Charge_Transfer
ilg_df.set_index('Conformer', inplace=True)
ilg_df
Out[479]:
In [498]:
il_pair_confs = pg.load_object('il_pair_analysis2')
In [500]:
emim_cl_df = il_pair_confs.get_table(filters={'Anion':'cl','Cation':'emim'}, mol=True)
emim_cl_df = emim_cl_df.set_index('Letter').loc[['2B', '2B', '2C', '2D', '2E']]
print emim_cl_df.columns
coro_cl_df = il_coro_confs.get_table(filters={'Cation':'emim', 'Anion':'cl'}, mol=True).sort(['Letter'])
coro_cl_df.set_index('Letter', inplace=True)
print coro_cl_df.columns
In [504]:
def diff(df1, df2, value1,value2=None):
if value2 is None: value2 = value1
return df1[value1]-df2[value2]
emim_hs = [(2,7), (4,6), (5,5), ('6a',17), ('6b',18), ('6c',19),
('7a',11), ('7b',12), ('8a',13), ('8b',14), ('8c',15)]
columns=['Energy', 'an_cat_dist', 'H-bond cat-an', 'C2_charge','Charge_Transfer']
for name, number in emim_hs:
label = 'H{}_charge'.format(name)
columns.append(label)
d_lst=[]
for (ix, df1), (ix2, df2) in zip(coro_cl_df.iterrows(), emim_cl_df.iterrows()):
d = [
diff(df1, df2, 'Energy of Association (kJmol^{-1}) Corrected'),
diff(df1, df2,'an_cat_dist' ,'Dist'),
diff(df1, df2,'H-bond cat-an','H-bond'),
diff(df1, df2,'C2_charge', 'C2_charge'),
diff(df1, df2,'Charge_Transfer')
]
for name, number in emim_hs:
label = 'H{}_charge'.format(name)
d.append(diff(df1, df2,label,label))
d_lst.append(d)
pd.DataFrame(d_lst, columns=columns).set_index(coro_cl_df.index)
Out[504]:
In [483]:
mol_path = os.path.join(inpath, 'Coronene-Ion')
coro_mol, cl_mol, bf_mol, etso_mol = [
pg.Molecule(os.path.join(inpath, 'Coronene-Ion'), init_fname='CJS1_{}_-_init.com'.format(m),
opt_fname='CJS1_{}_-_6-311+g-d-p-_gd3bj_opt_*.log'.format(m),
freq_fname='CJS1_{}_-_6-311+g-d-p-_gd3bj_freq_.log'.format(m),
nbo_fname='CJS1_{}_-_6-311+g-d-p-_-_pop-nbo-full-_*.log'.format(m))
for m in ['coro', 'cl', 'bf', 'etso']]
m='emim'
emim_mol = pg.Molecule(os.path.join(inpath, 'Coronene-Ion'), init_fname='CJS1_{}_-_init.com'.format(m),
opt_fname='CJS1_{}_-_6-311+g-d-p-_gd3bj_opt_*.log'.format(m),
freq_fname='CJS1_{}_-_6-311+g-d-p-_gd3bj_freq_.log'.format(m),
nbo_fname='CJS1_{}_-_6-311+g-d-p-_*_pop-nbo-full-_*.log'.format(m))
In [484]:
def plot_dos(mols, titles, ubound=10, lbound=-20, ufreq=10, linewidth=2, gap=True):
fig_dos, axes = plt.subplots(1,len(mols), sharex=True)
for mol, ax, title in zip(mols, axes, titles):
mol.plot_dos(ubound=ubound, lbound=lbound, ax=ax)
ax.set_title(title, {'fontsize': 8})
ax.tick_params(labelleft='off')
ax.set_xticks([0,2,4,6,8])
l=ax.legend()
l.set_visible(False)
for mol_ax, mol, mol_name in zip(axes, mols, titles):
homo, lumo = mol.get_orbital_energies(mol.get_orbital_homo_lumo())
if ufreq:
xlower, xupper = 0,ufreq
else:
xlower, xupper = mol_ax.get_xbound()
mol_ax.plot([xlower, xupper], [homo,homo], 'g-', linewidth=linewidth)
mol_ax.plot([xlower, xupper], [lumo,lumo], 'r-', linewidth=linewidth)
mol_ax.annotate('{}eV'.format(homo.round(1)), xy=(xupper, homo), xytext=(-5, -10), ha='right',
textcoords='offset points')
mol_ax.annotate('{}eV'.format(lumo.round(1)), xy=(xupper, lumo), xytext=(-5, 5), ha='right',
textcoords='offset points')
gap = lumo-homo
if gap:
mol_ax.annotate('{}eV'.format(gap.round(1)), xy=(xupper, homo+0.5*gap), xytext=(-5, -4), ha='right',
textcoords='offset points')
axes[0].set_ylabel('Energy (eV)')
axes[0].tick_params(labelleft='on')
if ufreq: axes[0].set_xbound(0,ufreq)
# The big subplot for x label
ax = fig_dos.add_subplot(111)
ax.tick_params(top='off', bottom='off', left='off', right='off',
labelbottom='on', labelleft='on', pad=25)
ax.set_xticklabels([])
ax.set_yticklabels([])
ax.set_frame_on(False)
ax.set_xlabel('Density of States (eV$^{-1}$)')
fig_dos.tight_layout(w_pad=0.1)
return fig_dos
fig_dos={}
for anion, anion_mol, anion_name in zip(['cl', 'bf', 'etso'], [cl_mol, bf_mol, etso_mol],
['[Cl]$^-$', '[BF$_4$]$^-$', '[C$_2$SO$_4$]$^-$']):
a,b=[coro_mol, emim_mol, anion_mol],['C$_{24}$H$_{12}$','[C$_2$C$_1$im]$^+$',anion_name]
for ix, df in il_coro_confs.get_table(filters={'Anion':anion,'Cation':'emim'}, mol=True).iterrows():
a.append(df.Molecule)
b.append(df.Letter)
fig_dos[anion] = plot_dos(a, b)
fig_dos[anion].set_size_inches(8,3)
#fig_dos[anion].savefig('test_'+anion, dpi=400)
#'C$_{24}$H$_{12}$','[C2C1im]$^+$','[Cl]$^-$'
In [485]:
il_coro_confs.add_mol_property_subset(['HOMO', 'LUMO'], 'get_orbital_homo_lumo')
h, l= coro_mol.get_orbital_energies(coro_mol.get_orbital_homo_lumo())
coro_bg = l-h
print coro_bg
for ix, df in il_coro_confs._df.iterrows():
h, l= df.Molecule.get_orbital_energies(df.Molecule.get_orbital_homo_lumo())
print df.Letter, l-h-coro_bg
In [486]:
h,l = il_coro_confs.get_molecule(23).get_orbital_homo_lumo()
print il_coro_confs.get_molecule(23).get_orbital_energies([l, h]).round(2)
h,l = il_coro_confs.get_molecule(24).get_orbital_homo_lumo()
print il_coro_confs.get_molecule(24).get_orbital_energies([l, h, h-6, h-28,h-42]).round(2)
h,l = il_coro_confs.get_molecule(46).get_orbital_homo_lumo()
print il_coro_confs.get_molecule(46).get_orbital_energies([l, h, h-11, h-14, h-54]).round(2)
In [487]:
doc = pg.MSDocument()
doc.add_heading('Ionic Liquid Summary', level=0)
doc.add_mpl(fig_ecl, dpi=180, width=12)
doc.add_mpl(fig_ebf, dpi=180, width=12)
doc.add_mpl(fig_eetso, dpi=180, width=12)
doc.add_dataframe(ilg_df, sig_figures=4)
doc.add_mpl(ax_energy.get_figure(), dpi=96, height=8)
doc.add_mpl(fig_dos['cl'], dpi=200, width=12)
doc.add_mpl(fig_dos['bf'], dpi=200, width=12)
doc.add_mpl(fig_dos['etso'], dpi=200, width=12)
doc.save('For_Final_Report_IL_Coro.docx')
In [526]:
mol = il_coro_confs.get_molecule(19)
mol.get_opt_energy(units='kJmol-1', zpe_correct=True)
Out[526]:
In [528]:
mol = il_coronene.get_molecule(9)
mol.get_opt_energy(units='kJmol-1', zpe_correct=True)
Out[528]: