Skip to content

Commit 017ab86

Browse files
tjduignatjduigna
committed
Computing energy curvature, disabled basis-related custom_traits
1 parent 99cf8b5 commit 017ab86

6 files changed

Lines changed: 188 additions & 144 deletions

File tree

‎exatomic/algorithms/delocalization.py‎

Lines changed: 174 additions & 137 deletions
Original file line numberDiff line numberDiff line change
@@ -1,150 +1,187 @@
11
# -*- coding: utf-8 -*-
22
# Copyright (c) 2015-2016, Exa Analytics Development Team
33
# Distributed under the terms of the Apache License 2.0
4-
'''
4+
"""
55
Delocalization
66
################################
7-
Miscellaneous functions for computing the delocalization error
8-
inherent in carefully constructed exatomic universes. This universe
9-
requires results from 3 different quantum chemical calculations on
7+
Miscellaneous functions for computing the curvature in the energy
8+
of a system as a function of the electron number. These functions
9+
require results from 3 different quantum chemical calculations on
1010
an (N-1), N, and (N+1) electron system.
11-
'''
11+
"""
1212

13+
import numpy as np
14+
import pandas as pd
1315

14-
def stack_unis(unis, same_frame=False):
15-
'''
16-
A non-general algorithm for combining "universes" into a single universe.
17-
'''
18-
attrs = [key[1:] for key, value in vars(unis[0]).items() if key not in
19-
['_traits_need_update', '_widget', '_lines'] and key.startswith('_')]
20-
dfs = {}
21-
if same_frame:
22-
for attr in attrs:
23-
if attr == 'field':
24-
continue
25-
dfs[attr] = getattr(unis[0], attr)
26-
if hasattr(unis[0], 'field'):
27-
field_values = [i.field.field_values[0] for i in unis]
28-
field_params = pd.concat([pd.DataFrame(unis[0].field)] * len(unis)).reset_index(drop=True)
29-
field_params['frame'] = 0
30-
#field_params['label'] = field_params['label'].astype(np.int64)
31-
dfs['field'] = exatomic.field.AtomicField(field_params, field_values=field_values)
32-
return exatomic.Universe(**dfs)
33-
for attr in attrs:
34-
if attr == 'field':
35-
continue
36-
subdfs = []
37-
for i, uni in enumerate(unis):
38-
smdf = getattr(uni, attr)
39-
smdf['frame'] = i
40-
subdfs.append(smdf)
41-
adf = pd.concat(subdfs)
42-
adf.reset_index(drop=True, inplace=True)
43-
dfs[attr] = adf
44-
if hasattr(unis[0], 'field'):
45-
field_values = [i.field.field_values[0] for i in unis]
46-
field_params = pd.concat([pd.DataFrame(unis[0].field)] * len(unis)).reset_index(drop=True)
47-
field_params['frame'] = range(len(unis))
48-
field_params['label'] = 0
49-
field_params['label'] = field_params['label'].astype(np.int64)
50-
dfs['field'] = exatomic.field.AtomicField(field_params, field_values=field_values)
51-
return exatomic.Universe(**dfs)
16+
from exatomic import Energy
17+
#from exatomic.field import AtomicField
18+
#from exatomic.container import Universe
5219

20+
#def stack_unis(unis, same_frame=False):
21+
# """
22+
# A non-general algorithm for combining ``universes'' into a single universe.
23+
# """
24+
# attrs = [key[1:] for key, value in vars(unis[0]).items() if key not in
25+
# ['_traits_need_update', '_widget', '_lines'] and key.startswith('_')]
26+
# dfs = {}
27+
# if same_frame:
28+
# for attr in attrs:
29+
# if attr == 'field':
30+
# continue
31+
# dfs[attr] = getattr(unis[0], attr)
32+
# if hasattr(unis[0], 'field'):
33+
# field_values = [i.field.field_values[0] for i in unis]
34+
# field_params = pd.concat([pd.DataFrame(unis[0].field)] * len(unis)).reset_index(drop=True)
35+
# field_params['frame'] = 0
36+
# #field_params['label'] = field_params['label'].astype(np.int64)
37+
# dfs['field'] = AtomicField(field_params, field_values=field_values)
38+
# return Universe(**dfs)
39+
# for attr in attrs:
40+
# if attr == 'field':
41+
# continue
42+
# subdfs = []
43+
# for i, uni in enumerate(unis):
44+
# smdf = getattr(uni, attr)
45+
# smdf['frame'] = i
46+
# subdfs.append(smdf)
47+
# adf = pd.concat(subdfs)
48+
# adf.reset_index(drop=True, inplace=True)
49+
# dfs[attr] = adf
50+
# if hasattr(unis[0], 'field'):
51+
# field_values = [i.field.field_values[0] for i in unis]
52+
# field_params = pd.concat([pd.DataFrame(unis[0].field)] * len(unis)).reset_index(drop=True)
53+
# field_params['frame'] = range(len(unis))
54+
# field_params['label'] = 0
55+
# field_params['label'] = field_params['label'].astype(np.int64)
56+
# dfs['field'] = AtomicField(field_params, field_values=field_values)
57+
# return Universe(**dfs)
5358

54-
class Deloc:
55-
@staticmethod
56-
def compute_deloc(uni, debug, jtype=None):
57-
#Basics
58-
cat = uni.frames[0]
59-
neut = uni.frames[1]
60-
an = uni.frames[2]
61-
cat_en = cat.energy
62-
neut_en = neut.energy
63-
an_en = an.energy
64-
#Electron counts
65-
aco = cat.orbitals.get_orbital()
66-
acate = cat.orbitals.get_orbital().orbital
67-
try:
68-
bco = cat.orbitals.get_orbital(spin=1)
69-
except IndexError:
70-
bco = aco
71-
ano = neut.orbitals.get_orbital()
72-
try:
73-
bno = neut.orbitals.get_orbital(spin=1)
74-
except IndexError:
75-
bno = ano
76-
aao = an.orbitals.get_orbital()
77-
try:
78-
bao = an.orbitals.get_orbital(spin=1)
79-
except IndexError:
80-
bao = aao
59+
def combine_deloc(delocs, order=None):
60+
"""
61+
Given a list of the results of compute_deloc, return a single
62+
dataframe containing all of the E(N) results.
63+
"""
64+
if order is not None:
65+
oldorder = [deloc.name for deloc in delocs]
66+
reordered = []
67+
for ordr in order:
68+
for i, old in enumerate(oldorder):
69+
if ordr.upper() == old.upper():
70+
reordered.append(delocs[i])
71+
else:
72+
reordered = delocs
73+
for i, deloc in enumerate(reordered):
74+
if i > 0:
75+
reordered[i].drop('n', axis=1, inplace=True)
76+
return pd.concat(reordered, axis=1)
8177

82-
#Eigenvalues
83-
if aco.orbital < ano.orbital:
84-
lumoca = cat.orbitals.get_orbital(index=aco.name + 1).energy
85-
homone = ano.energy
86-
if bco.orbital < bno.orbital:
87-
lumoca = cat.orbitals.get_orbital(index=bco.name + 1).energy
88-
homone = bno.energy
89-
if ano.orbital < aao.orbital:
90-
lumone = neut.orbitals.get_orbital(index=ano.name + 1).energy
91-
homoan = aao.energy
92-
if bno.orbital < bao.orbital:
93-
lumone = neut.orbitals.get_orbital(index=bno.name + 1).energy
94-
homoan = bao.energy
9578

96-
#Compute J^2
97-
jone = homone + (cat_en - neut_en)
98-
jtwo = homoan + (neut_en - an_en)
99-
jtype = None
100-
j2 = jone ** 2 + jtwo ** 2
101-
if jtype == 'EA':
102-
j2 = jone ** 2
103-
elif jtype == 'IP':
104-
j2 = jtwo ** 2
105-
#Compute E(n) and curvature coefficients
106-
q = np.linspace(0, 1, 51)
107-
negdE = an_en - neut_en
108-
posdE = neut_en - cat_en
109-
autoev = 27.2114
110-
negE = (negdE*q + ((lumone - negdE)*(1 - q) + (negdE - homoan)*q)*q*(1 - q)) * autoev
111-
posE = (-posdE*q + ((lumoca - posdE)*(1 - q) + (posdE - homone)*q)*q*(1 - q)) * autoev
112-
ancur = (np.sum((lumone - negdE)*(6*q - 4) + (negdE - homoan)*(2 - 6*q)))/(len(q)*2) * autoev
113-
catcur = (np.sum((lumoca - posdE)*(6*q - 4) + (posdE - homone)*(2 - 6*q)))/(len(q)*2) * autoev
114-
data = np.empty((len(q)*2 - 1,), dtype = [('n', 'f8'), ('E', 'f8')])
115-
data['n'] = np.concatenate((np.fliplr([-q])[0][:-1], q))
116-
data['E'] = np.concatenate((np.fliplr([posE])[0][:-1], negE))
117-
#Proper object and tack on tidbits
118-
retvar = pd.DataFrame(data)
119-
retvar.ancur = ancur
120-
retvar.catcur = catcur
121-
retvar.j2 = j2
122-
retvar.name = uni.name
123-
retvar.focc = uni.focc
124-
if debug:
125-
print('============', uni.name, '============')
126-
print('aco = \n', aco)
127-
print('bco = \n', bco)
128-
if aco.orbital < ano.orbital:
129-
print('lumoca = \n', cat.orbitals.get_orbital(index=aco.name + 1))
130-
if bco.orbital < bno.orbital:
131-
print('lumoca = \n', cat.orbitals.get_orbital(index=bco.name + 1))
132-
print('ano = \n', ano)
133-
print('bno = \n', bno)
134-
if ano.orbital < aao.orbital:
135-
print('lumone = \n', neut.orbitals.get_orbital(index=ano.name + 1))
136-
if bno.orbital < bao.orbital:
137-
print('lumone = \n', neut.orbitals.get_orbital(index=bno.name + 1))
138-
print('aao = \n', aao)
139-
print('bao = \n', bao)
140-
print('lumoca energy = ', lumoca)
141-
print('homone energy = ', homone)
142-
print('lumone energy = ', lumone)
143-
print('homoan energy = ', homoan)
144-
print('cat energy = ', cat_en)
145-
print('neut energy = ', neut_en)
146-
print('an energy = ', an_en)
147-
return retvar
79+
def compute_deloc(cat, neut, an, tag='', debug=False, jtype=None):
80+
"""
81+
Computes the curvature of the energy of a system as a function
82+
of the number of electrons in the system E(N).
14883
149-
def __new__(self, universe, debug=None, **kwargs):
150-
return self.compute_deloc(universe, debug=debug, **kwargs)
84+
Args
85+
cat (exatomic.Universe): N-1 electron system
86+
neut (exatomic.Universe): N electron system
87+
an (exatomic.Universe): N+1 electron system
88+
89+
Returns
90+
ret (pd.DataFrame): The energy as a function of N
91+
"""
92+
cat_en = cat.frame.ix[0].total_energy
93+
neut_en = neut.frame.ix[0].total_energy
94+
an_en = an.frame.ix[0].total_energy
95+
96+
# Get the highest occupied molecular orbitals of each system
97+
alpha_cat_homo = cat.orbital.get_orbital()
98+
alpha_neut_homo = neut.orbital.get_orbital()
99+
alpha_an_homo = an.orbital.get_orbital()
100+
101+
# Check for open shell nature of any of the systems
102+
try:
103+
beta_cat_homo = cat.orbital.get_orbital(spin=1)
104+
except IndexError:
105+
beta_cat_homo = alpha_cat_homo
106+
try:
107+
beta_neut_homo = neut.orbital.get_orbital(spin=1)
108+
except IndexError:
109+
beta_neut_homo = alpha_neut_homo
110+
try:
111+
beta_an_homo = an.orbital.get_orbital(spin=1)
112+
except IndexError:
113+
beta_an_homo = alpha_an_homo
114+
115+
# Find the right orbital energies
116+
if alpha_cat_homo.vector < alpha_neut_homo.vector:
117+
lumoca = cat.orbital.get_orbital(index=alpha_cat_homo.name + 1).energy
118+
homone = alpha_neut_homo.energy
119+
if beta_cat_homo.vector < beta_neut_homo.vector:
120+
lumoca = cat.orbital.get_orbital(index=beta_cat_homo.name + 1).energy
121+
homone = beta_neut_homo.energy
122+
if alpha_neut_homo.vector < alpha_an_homo.vector:
123+
lumone = neut.orbital.get_orbital(index=alpha_neut_homo.name + 1).energy
124+
homoan = alpha_an_homo.energy
125+
if beta_neut_homo.vector < beta_an_homo.vector:
126+
lumone = neut.orbital.get_orbital(index=beta_neut_homo.name + 1).energy
127+
homoan = beta_an_homo.energy
128+
129+
#Compute J^2
130+
jone = homone + (cat_en - neut_en)
131+
jtwo = homoan + (neut_en - an_en)
132+
jtype = None
133+
j2 = jone ** 2 + jtwo ** 2
134+
if jtype == 'EA':
135+
j2 = jone ** 2
136+
elif jtype == 'IP':
137+
j2 = jtwo ** 2
138+
139+
#Compute E(n) and curvature coefficients
140+
q = np.linspace(0, 1, 51)
141+
negdE = an_en - neut_en
142+
posdE = neut_en - cat_en
143+
autoev = Energy['Ha', 'eV']
144+
negE = (negdE*q + ((lumone - negdE)*(1 - q) + (negdE - homoan)*q)*q*(1 - q)) * autoev
145+
posE = (-posdE*q + ((lumoca - posdE)*(1 - q) + (posdE - homone)*q)*q*(1 - q)) * autoev
146+
ancur = (np.sum((lumone - negdE)*(6*q - 4) + (negdE - homoan)*(2 - 6*q)))/(len(q)*2) * autoev
147+
catcur = (np.sum((lumoca - posdE)*(6*q - 4) + (posdE - homone)*(2 - 6*q)))/(len(q)*2) * autoev
148+
colname = '{} ({:.2f},{:.2f})'.format(tag, catcur, ancur)
149+
data = np.empty((len(q)*2 - 1,), dtype = [('n', 'f8'), (colname, 'f8')])
150+
data['n'] = np.concatenate((np.fliplr([-q])[0][:-1], q))
151+
data[colname] = np.concatenate((np.fliplr([posE])[0][:-1], negE))
152+
153+
#Proper object and tack on tidbits
154+
ret = pd.DataFrame(data)
155+
ret.ancur = ancur
156+
ret.catcur = catcur
157+
ret.j2 = j2
158+
ret.name = tag
159+
ret.colname = colname
160+
if debug:
161+
print('============', tag, '============')
162+
print('alpha cation HOMO =', alpha_cat_homo, sep='\n')
163+
print('beta cation HOMO =', beta_cat_homo, sep='\n')
164+
if alpha_cat_homo.vector < alpha_neut_homo.vector:
165+
print('lumoca =',
166+
cat.orbital.get_orbital(index=alpha_cat_homo.name + 1), sep='\n')
167+
if beta_cat_homo.vector < beta_neut_homo.vector:
168+
print('lumoca =',
169+
cat.orbital.get_orbital(index=beta_cat_homo.name + 1), sep='\n')
170+
print('alpha neutral HOMO =', alpha_neut_homo, sep='\n')
171+
print('beta neutral HOMO =', beta_neut_homo, sep='\n')
172+
if alpha_neut_homo.vector < alpha_an_homo.vector:
173+
print('lumone =',
174+
neut.orbital.get_orbital(index=alpha_an_homo.name + 1), sep='\n')
175+
if beta_neut_homo.vector < beta_an_homo.vector:
176+
print('lumone =',
177+
neut.orbital.get_orbital(index=beta_an_homo.name + 1), sep='\n')
178+
print('alpha anion HOMO =', alpha_an_homo, sep='\n')
179+
print('beta anion HOMO =', beta_an_homo, sep='\n')
180+
print('lumoca energy = ', lumoca)
181+
print('homone energy = ', homone)
182+
print('lumone energy = ', lumone)
183+
print('homoan energy = ', homoan)
184+
print('cat energy = ', cat_en)
185+
print('neut energy = ', neut_en)
186+
print('an energy = ', an_en)
187+
return ret

‎exatomic/algorithms/geometry.py‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -13,10 +13,10 @@
1313

1414
def make_small_molecule(center=None, ligand=None, distance=None, geometry=None,
1515
offset=None, plane=None, axis=None, domains=None, unit='A'):
16-
'''
16+
"""
1717
A minimal molecule builder for simple one-center, homogeneous ligand
18-
molecules of various general chemistry molecular geometries. If domains
19-
is not specified and geometry is ambiguous (like 'bent'),
18+
molecules of various general chemistry molecular geometries. If `domains'
19+
is not specified and geometry is ambiguous (like `bent'),
2020
it just guesses the simplest geometry (smallest number of domains).
2121
2222
Args
@@ -31,7 +31,7 @@ def make_small_molecule(center=None, ligand=None, distance=None, geometry=None,
3131
3232
Returns
3333
exatomic.atom.Atom: Atom table of small molecule
34-
'''
34+
"""
3535
if center is None or ligand is None or distance is None or geometry is None:
3636
raise NotImplementedError('must supply center, ligand, distance and geometry')
3737
distance *= Length[unit, 'au']

‎exatomic/atom.py‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -15,6 +15,7 @@
1515
symbol_to_element_mass)
1616
from exatomic.error import PeriodicUniverseError
1717
from exatomic.algorithms.distance import minimal_image_counts
18+
from exatomic.algorithms.geometry import make_small_molecule
1819

1920

2021
class Atom(DataFrame):

‎exatomic/basis.py‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -135,7 +135,7 @@ class GaussianBasisSet(BasisSet):
135135
_index = 'primitive'
136136
_traits = ['shell_function']
137137
_precision = {'alpha': 8, 'd': 8}
138-
_categories = {'set': np.int64, 'l': np.int64, 'shell_function': np.int64}
138+
_categories = {'set': np.int64, 'L': np.int64, 'shell_function': np.int64}
139139

140140
# def _custom_traits(self):
141141
# g = self.grpd

‎exatomic/container.py‎

Lines changed: 0 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -41,7 +41,6 @@ class Meta(TypedMeta):
4141
frame = Frame
4242
atom_two = AtomTwo
4343
unit_atom = UnitAtom
44-
periodic_two = PeriodicTwo
4544
projected_atom = ProjectedAtom
4645
visual_atom = VisualAtom
4746
molecule = Molecule

‎exatomic/orbital.py‎

Lines changed: 8 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -56,11 +56,18 @@ class Orbital(DataFrame):
5656
Note:
5757
Spin zero means alpha spin or unknown and spin one means beta spin.
5858
"""
59-
_columns = ['frame', 'energy', 'x', 'y', 'z', 'occupation', 'spin', 'vector']
59+
_columns = ['frame', 'energy', 'occupation', 'spin', 'vector']
6060
_index = 'orbital'
6161
_groupby = ('frame', np.int64)
6262
_categories = {'spin': np.int64}
6363

64+
def get_orbital(self, orb=-1, spin=0, index=None):
65+
if index is None:
66+
return self[(self['occupation'] > 0) & (self['spin'] == spin)].iloc[orb]
67+
else:
68+
return self.iloc[index]
69+
70+
6471

6572
class MOMatrix(DataFrame):
6673
"""

0 commit comments

Comments
 (0)