Skip to content

Commit

Permalink
cleaning up splien fitter
Browse files Browse the repository at this point in the history
  • Loading branch information
@snelliott
snelliott committed Oct 16, 2021
1 parent 6e544b8 commit dfb6b64
Show file tree
Hide file tree
Showing 2 changed files with 101 additions and 58 deletions.
152 changes: 98 additions & 54 deletions automol/pot/_fit.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -5,14 +5,12 @@

import numpy
from scipy.interpolate import interp1d
from scipy.interpolate import CubicSpline


def fit_1d_potential(pot_dct, min_thresh=-0.0001, max_thresh=50.0):
""" Get a physical hindered rotor potential via a series of spline fits
"""

pot = list(pot_dct.values())

def _initialize_pot(pot_dct):
grid_vals, pot = pot_dct.items()
pot = list(pot)
# Check if all values are bad
if all(numpy.isclose(-10.0, val) for val in pot):
print('ERROR: All potential values missing')
Expand All @@ -21,7 +19,15 @@ def fit_1d_potential(pot_dct, min_thresh=-0.0001, max_thresh=50.0):
lpot = len(pot)+1
pot.append(0.0)

return grid_vals, pot, lpot


def _round_near_zero(pot, max_thresh, min_thresh):
""" Caps the potentials over a given threshold
"""
# Print warning messages
# Build a potential list from only successful calculations
# First replace high potential values with max_thresh
print_pot = False
if any(val > max_thresh for val in pot):
print_pot = True
Expand All @@ -32,7 +38,9 @@ def fit_1d_potential(pot_dct, min_thresh=-0.0001, max_thresh=50.0):
# reset any negative values for the first grid point to 0.
if pot[0] < 0.:
print('ERROR: The first potential value should be 0.')
pot[0] = 0.
pot[0] = 0.0
elif pot[0] < .01:
pot[0] = 0.0
if any(val < min_thresh for val in pot):
print_pot = True
print(f'Warning: Found pot val of {min(pot):.2f}',
Expand All @@ -41,86 +49,122 @@ def fit_1d_potential(pot_dct, min_thresh=-0.0001, max_thresh=50.0):

if print_pot:
print('Potential before spline:', pot)
return pot

# Build a potential list from only successful calculations
# First replace high potential values with max_thresh
# Then replace any negative potential values cubic spline fit values
idx_success = []
pot_success = []

def _cap_max_thresh_drop_neg(pot, lpot, min_thresh, max_thresh):
x_idxs = []
y_pots = []
for idx in range(lpot):
if pot[idx] < 600. and pot[idx] > min_thresh:
idx_success.append(idx)
x_idxs.append(idx)
if pot[idx] < max_thresh:
pot_success.append(pot[idx])
y_pots.append(pot[idx])
else:
pot_success.append(max_thresh)

if len(pot_success) > 3:
# Build a new potential list using a spline fit of the HR potential

spline_vals = pot_success[:-1]
mid_idx = round(len(spline_vals)/2)
spline_vals = spline_vals[mid_idx:] + spline_vals[:mid_idx]
pot_spl = interp1d(
numpy.array(idx_success[:-1]), numpy.array(spline_vals), kind='cubic')
for idx in range(lpot):
pot[idx] = float(pot_spl(idx))
if not mid_idx % 2:
pot = pot[mid_idx+1:] + pot[:mid_idx+1]
else:
pot = pot[mid_idx:] + pot[:mid_idx]
pot.append(0.0)
y_pots.append(max_thresh)
return numpy.array(x_idxs), numpy.array(y_pots)


def _periodic_spline(x_vals, y_vals, orig_len):
""" perform a cubic periodic spline on an x and y array
"""
spl_fit = []
if len(y_vals) > 3:
spl_fun = CubicSpline(x_vals, y_vals, bc_type='periodic')
for idx in range(orig_len):
spl_fit.append(float(spl_fun(idx)))
return spl_fit


def _spline_cap_max_thresh(pot, lpot, max_thresh, min_thresh):
"""replace any negative potential values cubic spline fit values
"""
x_idxs, y_pots = _cap_max_thresh_drop_neg(
pot, lpot, min_thresh, max_thresh)
spl_fit = _periodic_spline(x_idxs, y_pots, lpot)
return pot if spl_fit else spl_fit


def _spline_no_alteration(pot, lpot, min_thresh):
"""Do a spline fit for negatives without changing anything else
"""
x_pos = numpy.array([i for i in range(lpot)
if pot[i] >= min_thresh])
y_pos = numpy.array([pot[i] for i in range(lpot)
if pot[i] >= min_thresh])
spl_fit = _periodic_spline(x_pos, y_pos, lpot)
return pot if spl_fit else spl_fit


def fit_1d_potential(pot_dct, min_thresh=-0.0001, max_thresh=50.0):
""" Get a physical hindered rotor potential via a series of spline fits
"""
# unpack potential dictionary
grid_vals, pot, lpot = _initialize_pot(pot_dct)
pot = _round_near_zero(pot, max_thresh, min_thresh)
# Build a new potential list using a spline fit of the HR potential
pot = _spline_cap_max_thresh(pot, lpot, max_thresh, min_thresh)

# Do second spline fit of only positive values if any negative values found
if any(val < min_thresh for val in pot):
print('Still found negative potential values after first spline')
print('Potential after spline:', pot)
if len(pot_success) > 3:
x_pos = numpy.array([i for i in range(lpot)
if pot[i] >= min_thresh])
y_pos = numpy.array([pot[i] for i in range(lpot)
if pot[i] >= min_thresh])
pos_pot_spl = interp1d(x_pos, y_pos, kind='cubic')
pot_pos_fit = []
for idx in range(lpot):
pot_pos_fit.append(pos_pot_spl(idx))
else:
pot_pos_fit = []
for idx in range(lpot):
pot_pos_fit.append(pot[idx])

print('Potential after spline:', pot_pos_fit)
print('Potential after first spline:', pot)
pot = _spline_no_alteration(pot, lpot, min_thresh)
print('Potential after second spline:', pot)
# Perform second check to see if negative potentials have been fixed
if any(val < min_thresh for val in pot_pos_fit):
if any(val < min_thresh for val in pot):
print('Still found negative potential values after second spline')
print('Replace with linear interpolation of positive values')
neg_idxs = [i for i in range(lpot) if pot_pos_fit[i] < min_thresh]
neg_idxs = [i for i in range(lpot) if pot[i] < min_thresh]
clean_pot = []
for i in range(lpot):
if i in neg_idxs:
# Find the indices for positive vals around negative value
idx_0 = i - 1
idx_1 = None
while idx_0 in neg_idxs:
idx_0 = idx_0 - 1
for j in range(i, lpot):
if pot_pos_fit[j] >= min_thresh:
if pot[j] >= min_thresh:
idx_1 = j
break
if idx_1 is None:
for j in range(0, idx_0):
idx_1 = j + lpot
# Get a new value for this point on the potential by
# doing a linear interp of positives
interp_val = (
pot_pos_fit[idx_0] * (1.0-((i-idx_0)/(idx_1-idx_0))) +
pot_pos_fit[idx_1] * ((i-idx_0)/(idx_1-idx_0))
pot[idx_0] * (1.0-((i-idx_0)/(idx_1-idx_0))) +
(pot[idx_1 if idx_1 < lpot else idx_1 - lpot]
* ((i-idx_0)/(idx_1-idx_0)))
)
clean_pot.append(interp_val)
else:
clean_pot.append(pot_pos_fit[i])
clean_pot.append(pot[i])
final_potential = clean_pot.copy()
print('Potential after linear interp:', final_potential)

else:
final_potential = pot_pos_fit.copy()
final_potential = pot.copy()

# check for quadratic behavior around minimum
else:
final_potential = pot.copy()
#if lpot > 7:
# pos = final_potential[1]
# neg = final_potential[-2]
# if abs(pos - neg)/(pos + neg) > 0.3:
# print('Refitting potentials near the minimum to a quadratic')
# pos_1 = final_potential[2]
# neg_1 = final_potential[-3]
# quad_pots = [neg_1, neg, 0, pos, pos_1]
# quad_idxs = [0, 1, 2, 3, 4]
# quad_pot_coeffs = numpy.polyfit(quad_idxs, quad_pots, 2)
# print(quad_pots)
# print(quad_pot_coeffs)
# final_potential[-2] = quad_pot_coeffs[2] + quad_pot_coeffs[1] * 1 + quad_pot_coeffs[0] * 1**2
# final_potential[1] = quad_pot_coeffs[2] + quad_pot_coeffs[1] * 3 + quad_pot_coeffs[0] * 3**2
# print('Potential to fit quadratic', final_potential)

final_potential = final_potential[:-1]

Expand Down
7 changes: 3 additions & 4 deletions phydat/phycon.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -15,10 +15,9 @@
SOL = (qcc.get('speed of light in vacuum') *
qcc.conversion_factor('meter / second', 'bohr hartree / h'))
SOLMS = qcc.get('speed of light in vacuum')
KB = qcc.get('kb') # The Boltzmann constant (JK$^{-1}$)
H = qcc.get('h') # The Planck constant (Js)
HBAR = qcc.get('hbar') # The Planck constant (Js) over 2pi

KB = qcc.get('kb') # The Boltzmann constant (JK$^{-1}$)
H = qcc.get('h') # The Planck constant (Js)
HBAR = qcc.get('hbar') # The Planck constant (Js) over 2pi

# Energy Conversion factors
KCAL2CAL = qcc.conversion_factor('kcal/mol', 'cal/mol')
Expand Down

0 comments on commit dfb6b64

Please sign in to comment.