from __future__ import print_function
import numpy as np
import pandas as pd
from six import string_types
from ase import units
[docs]def expandrange(modestr):
'''
Convert a comma separated string of indices and dash separated ranges into
a list of integer indices
Parameters
----------
modestr : str
Returns
-------
indices : list of ints
Examples
--------
>>> from panther.pes import expandrange
>>> s = "2,3,5-10,20,25-30"
>>> expandrange(s)
[2, 3, 5, 6, 7, 8, 9, 10, 20, 25, 26, 27, 28, 29, 30]
'''
ranges = (x.split("-") for x in modestr.split(","))
return [i for r in ranges for i in range(int(r[0]), int(r[-1]) + 1)]
[docs]def calculate_energies(images, calc, modes='all'):
'''
Given a set of images as a nested OrderedDict of Atoms objects and a
calculator, calculate the energy for each displaced structure
Parameters
----------
images : OrderedDict
A nested OrderedDict of displaced Atoms objects
calc : calculator instance
ASE calculator
modes : str or list
Mode for which the PES will be calculated
Returns
-------
energies : pandas.DataFrame
DataFrame with the energies per displacement
'''
if isinstance(modes, string_types):
if modes.lower() in ['all', ':']:
modes = range(len(images))
else:
modes = expandrange(modes)
elif isinstance(modes, (list, tuple)):
pass
else:
ValueError('<modes> should be a str, list or tuple '
'got: {}'.format(type('modes')))
ecols = ['E_' + str(i) for i in range(-4, 5)]
energies = pd.DataFrame(0.0, columns=ecols, index=pd.Index(modes,
name='mode'))
for mode in modes:
for point in images[mode].keys():
print('# calculating energy for mode: {} point: {}'.format(mode,
point))
atoms = images[mode][point]
atoms.set_calculator(calc)
energy = atoms.get_potential_energy()
print('E[{0:d}, {1:d}] : {2:25.12f}'.format(mode, point, energy))
energies.set_value(mode, 'E_' + str(point), energy)
return energies
[docs]def fit_potentials(modeinfo, energies):
'''
Fit the potentials with 6th and 4th order polynomials
Parameters
----------
modeinfo : pandas.DataFrame
DataFrame with per mode characteristics, displacements, masses
and a flag to mark it a mode is a stretching mode or not
energies : pd.DataFrame
Energies per displacement
Returns
-------
out : (coeffs6o, coeffs4o)
DataFrames with 6th and 4th polynomial coefficients fitted to the
potential
'''
energies = energies.subtract(energies['E_0'], axis=0)
energies = energies / units.Hartree
E = energies.as_matrix()
D = np.dot(modeinfo['displacement'].values.reshape(-1, 1),
np.arange(-4, 5).reshape(1, -1))
D = D / np.sqrt(modeinfo['effective_mass'].values).reshape(-1, 1)
D = D.astype(float)
coeffs6o = pd.DataFrame(index=energies.index,
columns=['c_' + str(x)
for x in np.arange(6, -1, -1)])
coeffs4o = pd.DataFrame(index=energies.index,
columns=['c_' + str(x)
for x in np.arange(4, -1, -1)])
for i in energies.index:
coeffs4o.loc[i] = np.polyfit(D[i], E[i], deg=4)
coeffs6o.loc[i] = np.polyfit(D[i], E[i], deg=6)
return coeffs6o, coeffs4o
[docs]def differentiate(displacements, energies, order=2):
'''
Calculate numerical detivatives using the central difference formula
Parameters
----------
displacements : array_like
energies : DataFrame
order : int
Order of the derivative
Notes
-----
Central difference coefficients taken from [1]_
.. [1] Fornberg, B. (1988). Generation of finite difference
formulas on arbitrarily spaced grids. Mathematics of Computation,
51(184), 699-699. doi:10.1090/S0025-5718-1988-0935077-0
'''
energies = energies.subtract(energies['E_0'], axis=0)
energies = energies / units.Hartree
E = energies.as_matrix()
cols = ['p2', 'p4', 'p6', 'p8']
dt = [('idx', int)] + list(zip(cols, [float] * 4))
coeffs1d = np.zeros(9, dtype=dt)
coeffs1d['idx'] = range(-4, 5)
coeffs1d['p2'] = [0.0, 0.0, 0.0, -0.5, 0.0, 0.5, 0.0, 0.0, 0.0]
coeffs1d['p4'] = [0.0, 0.0, 1.0/12.0, -2.0/3.0, 0.0, 2.0/3.0, -1.0/12.0, 0.0, 0.0]
coeffs1d['p6'] = [0.0, -1.0/60.0, 3.0/20.0 , -3.0/4.0, 0.0, 3.0/4.0, -3.0/20.0, 1.0/60.0, 0.0]
coeffs1d['p8'] = [1.0/280.0, -4.0/105.0, 1.0/5.0, -4.0/5.0, 0.0, 4.0/5.0, -1.0/5.0, 4.0/105.0, -1.0/280.0]
coeffs2d = np.zeros(9, dtype=dt)
coeffs2d['idx'] = range(-4, 5)
coeffs2d['p2'] = [0.0, 0.0, 0.0, 1.0, -2.0, 1.0, 0.0, 0.0, 0.0]
coeffs2d['p4'] = [0.0, 0.0, -1.0/12.0, 4.0/3.0, -5.0/2.0, 4.0/3.0, -1.0/12.0, 0.0, 0.0]
coeffs2d['p6'] = [0.0, 1.0/90.0, -3.0/20.0, 3.0/2.0, -49.0/18.0, 3.0/2.0, -3.0/20.0, 1.0/90.0, 0.0]
coeffs2d['p8'] = [-1.0/560.0, 8.0/315.0, -1.0/5.0, 8.0/5.0, -205.0/72.0, 8.0/5.0, -1.0/5.0, 8.0/315.0, -1.0/560.0]
if order == 1:
C = coeffs1d[cols].view(np.float64).reshape(coeffs1d.shape + (-1,))
return np.dot(E, C) / displacements[:, np.newaxis]
elif order == 2:
C = coeffs2d[cols].view(np.float64).reshape(coeffs2d.shape + (-1,))
return np.dot(E, C) / np.power(displacements, 2)[:, np.newaxis]
else:
raise NotImplementedError('{} order derivatives not available'.format(order))
[docs]def harmonic_potential(x, freq, mu):
'''
Calculate the harmonic potential
Parameters
----------
x : float of numpy.array
Coordinate
mu : float
Reduced mass
freq : float
Frequency in cm^-1
'''
# conversion factor from cm^-1 to hartrees
conv = 100.0 * units.J * units._hplanck * units._c / units.Hartree
kconst = mu * (freq * conv)**2
return 0.5 * kconst * x**2