Fitting_with_scipyandastropy.py

def OPLS_pot(R,SIG=10,EPS=0.05):
    sij=np.sqrt(SIG*3.15061)
    eij=np.sqrt(EPS*0.1521)
    return 4*eij*((sij/R)**12-(sij/R)**6)
xdata = Ri
ydata = pot(Ri)
popt, pcov = curve_fit(OPLS_pot, xdata, ydata)
print popt


######
import numpy as np
from astropy.modeling.models import custom_model
from astropy.modeling.fitting import LevMarLSQFitter

# Define model
@custom_model
def sum_of_gaussians(x, amplitude1=1., mean1=-1., sigma1=1.,
                        amplitude2=1., mean2=1., sigma2=1.):
    return (amplitude1 * np.exp(-0.5 * ((x - mean1) / sigma1)**2) +
            amplitude2 * np.exp(-0.5 * ((x - mean2) / sigma2)**2))
@custom_model
def LJ_pot(R,AA=10,BB=1.5):
    return ((AA/R**12)+(BB/R**6))
# Generate fake data
np.random.seed(0)
x = Ri
y = pot(Ri)

# Fit model to data
#m_init = sum_of_gaussians()
m_init = LJ_pot()
fit = LevMarLSQFitter()
m = fit(m_init, x, y)

# Plot the data and the best fit
plt.plot(x, y, 'o', color='k')
plt.plot(x, m(x), color='r', lw=2)
print fit.fit_info