2

I am trying to create the curve fit with scipy for the energy eigenvalues calculated from a 4x4 Hamiltonian matrix. In the following error "energies" corresponds to the function in which I define the Hamiltonian, "xdata" is an array given after and out of the function and corresponds to k and "e" is the energy eigenvalues that a get. The error seems to be at the Hamiltonian matrix. However if I run the code without the curve_fit everything works fine. I have also tried using np.array according to other questions I found here but again it doesn't work. If a give a specific xdata in the curve fit, like xdata[0], the code works but it doesn't help me much since I want the fit using all values. Does anyone know what is the problem? Thank you all in advance!


Traceback (most recent call last):
     File "fitest.py", line 70, in <module>
     popt, pcov = curve_fit(energies,xdata, e)#,
  File "/eb/software/Python/2.7.12-intel-2016b/lib/python2.7/site-packages/scipy/optimize/minpack.py", line 651, in curve_fit
    res = leastsq(func, p0, args=args, full_output=1, **kwargs)
  File "/eb/software/Python/2.7.12-intel-2016b/lib/python2.7/site-packages/scipy/optimize/minpack.py", line 377, in leastsq
    shape, dtype = _check_func('leastsq', 'func', func, x0, args, n)
  File "/eb/software/Python/2.7.12-intel-2016b/lib/python2.7/site-packages/scipy/optimize/minpack.py", line 26, in _check_func
    res = atleast_1d(thefunc(*((x0[:numinputs],) + args)))
  File "/eb/software/Python/2.7.12-intel-2016b/lib/python2.7/site-packages/scipy/optimize/minpack.py", line 453, in _general_function
    return function(xdata, *params) - ydata
  File "fitest.py", line 23, in energies
    [         0.0,              0.0,                      0.0,     ep-2*V4*cos(kpt*d) ]],dtype=complex)
TypeError: only length-1 arrays can be converted to Python scalars

Code:

from numpy import sin, cos, array

from scipy.optimize import curve_fit

from numpy import *

from numpy.linalg import *


def energies(kpt, a=1.0, b=2.0, c=3.0, f=4.0):

   e1=-15.0
   e2=-10.0

   d=1.0
   v0=(-2.0/d**2) 
   V1=a*v0
   V2=b*v0
   V3=c*v0
   V4=d*v0

   basis=('|S, s>', '|S,px>', '|S, py>', '|S,pz>')
   h=array([[ e1-2*V1*cos(kpt*d), -2*V2*1j*sin(kpt*d),   0.0,                   0.0 ], 
     [    2*V2*1j*sin(kpt*d), e2-2*V3*cos(kpt*d),       0.0,                        0.0],  
    [         0.0,              0.0,           e2-2*V4*cos(kpt*d),                 0.0],            
   [         0.0,              0.0,                      0.0,     e2-2*V4*cos(kpt*d) ]],dtype=complex) 

   e,psi=eigh(h)


   return e

print energies(kpt=0.0)
k2=0.4*2*pi/2.05
print energies(kpt=k2)


xdata = array([0.0,k2])
print xdata

popt, pcov = curve_fit(energies, xdata, e)


print " "
print popt
print " "
touki
  • 21
  • 5
  • 1
    How are we supposed to know, what is wrong with your code, when you don't show it? You should provide a [Minimal, Complete, and Verifiable example](https://stackoverflow.com/help/mcve). The most likely reason for this error message would be that your fit function uses the [math module instead of a numpy functions](https://stackoverflow.com/a/21697288/8881141) – Mr. T Feb 07 '18 at 12:40
  • Excuse me but this is the first time that I post something and I didn't know if the information given would suffice.I have added the code. The values here are random – touki Feb 07 '18 at 13:10

1 Answers1

1

Your problem has nothing to do with your fit, you run into the same problem, if you perform

print energies(xdata)

The reason for this error message is that you put an array kpt into h as an array element and then tell numpy, to transform this array kpt into a complex number. Numpy is kind enough to transform an array of length one into a scalar, which then can be transformed into a complex number. This explains, why you didn't get an error message with xdata[0]. You can easily reproduce your problem like this

import numpy as np
#all fine with an array of length one
xa = np.asarray([1])
a = np.asarray([[xa, 2, 3], [4, 5, 6]])
print a
print a.astype(complex)

#can't apply dtype = complex to an array with two elements
xb = np.asarray([1, 2])
b = np.asarray([[xb, 2, 3], [4, 5, 6]])
print b
print b.astype(complex)

Idk, what you were trying to achieve with your energies function, so I can only speculate, what you were aiming at, when constructing the h array. Maybe a 3D array like this?

kpt = np.asarray([1, 2, 3])

h = np.zeros(16 * len(kpt), dtype = complex).reshape(len(kpt), 4, 4)
h[:, 0, 0] = 2 * kpt + 1
h[:, 0, 1] = kpt ** 2
h[:, 3, 2] = np.sin(kpt)
print h
Mr. T
  • 11,960
  • 10
  • 32
  • 54
  • Dear Piinthesky, thank you for your answer.I tried what you proposed but it didn't fix the problem. So I looked it up a bit more and the code worked when I used the numpy.linalg.norm for the kpt. – touki Feb 12 '18 at 08:40
  • Good to hear. If you found a solution to your problem, it is encouraged to [answer your own question](https://stackoverflow.com/help/self-answer) and [accept it](https://stackoverflow.com/help/accepted-answer). – Mr. T Feb 12 '18 at 08:47