Chimie Théorique » QMX

# -*- coding: utf-8 -*-

from math import sqrt

import numpy as np

from ase.units import kJ

from ase.utils.eos import EquationOfState

def BulkModulus(images):

    v = np.array([abs(np.linalg.det(A)) for A in images.A])

    #import matplotlib.pyplot as plt

    import pylab as plt

    plt.ion()

    EquationOfState(v, images.E).plot()

"""

    fit = np.poly1d(np.polyfit(v**-(1.0 / 3), images.E, 3))

    fit1 = np.polyder(fit, 1)

    fit2 = np.polyder(fit1, 1)

    for t in np.roots(fit1):

        if t > 0 and fit2(t) > 0:

            break

    v0 = t**-3

    e0 = fit(t)

    B = t**5 * fit2(t) / 9 / kJ * 1.0e24  # Gpa

    import pylab

    import matplotlib

    #matplotlib.use('GTK')

    pylab.ion()

    x = 3.95

    pylab.figure(figsize=(x * 2.5**0.5, x))

    pylab.plot(v, images.E, 'o')

    x = np.linspace(min(v), max(v), 100)

    pylab.plot(x, fit(x**-(1.0 / 3)), '-r')

    pylab.xlabel(u'volume [Å^3]')

    pylab.ylabel(u'energy [eV]')

    pylab.title(u'E: %.3f eV, V: %.3f Å^3, B: %.3f GPa' % (e0, v0, B))

    pylab.show()

"""