RY_KBAR          = 10.0 * AU_GPA / 2.0 RYDBERG_SI       = HARTREE_SI/2.0   # J NA               = 6.02214129e23  # Avogadro constant, mol-1 def read_Etot(fname, ibrav=4, bc_as_a_ratio = True):
    """
    Read cell parameters *(a,b,c)* and the corresponding energies from input file *fname*. 
    Each set of cell parameters is stored in a numpy array of lenght 6 
    for *(a,b,c,alpha,beta,gamma)* respectively. This is done for a future possible 
    extension but for now only the first 3 elements are used (the others are always 0).
    All sets are stored in *celldmsx* and *Ex*, the former is a nE*6 matrix, 
    the latter is a nE array.
    
    *ibrav* identifies the Bravais lattice as in Quantum Espresso and is needed 
    in input (default is 4, i.e. hexagonal cell). The input file format depends
    on *ibrav*, for example in the hex case, the first two columns are for *a* and
    *c* and the third is for the energies.
    
    If *bc_as_a_ratio=True*, the input data are assumed to be given as 
    :math:`(a,b/a,c/a)` in the input file and hence converted into :math:`(a,b,c)`
    which is how they are always stored internally in :py:mod:`pyqha`.
    
    Units must be :math:`a.u.` and :math:`Ryd/cell`
 
    """
    celldmsx=[]
    Ex=[]
    
    if ibrav in (1,2,3):                # cubic systems (a,a,a)
        with open(fname, "r") as lines:
            for line in lines:
                celldms=np.zeros([6])
                linesplit=line.split()
                celldms[0]=float(linesplit[0])
                E=float(linesplit[1])
                celldmsx.append(celldms)
                Ex.append(E)
    
    elif ibrav in (4,6,7):              # hexagonal or tetragonal systems (a,a,c)
        with open(fname, "r") as lines:
            for line in lines:
                celldms=np.zeros([6])
                linesplit=line.split()
                celldms[0]=float(linesplit[0])
                if (bc_as_a_ratio==True):
                    celldms[2]=float(linesplit[1])*celldms[0]   # convert from c/a to c
                else:
                    celldms[2]=float(linesplit[1])              # as it is
                E=float(linesplit[2])
                celldmsx.append(celldms)
                Ex.append(E)
                
    elif ibrav in (8,9,10,11):          # orthorombic systems (a,b,c)
        with open(fname, "r") as lines:
            for line in lines:
                celldms=np.zeros([6])
                linesplit=line.split()
                celldms[0]=float(linesplit[0])
                if (bc_as_a_ratio==True):
                    celldms[1]=float(linesplit[1])*celldms[0]   # convert from b/a to b
                    celldms[2]=float(linesplit[2])*celldms[0]   # convert from c/a to c
                else:
                    celldms[1]=float(linesplit[1])              # as it is
                    celldms[2]=float(linesplit[2])              # as it is
                E=float(linesplit[3])
                celldmsx.append(celldms)
                Ex.append(E)    
    else:
        return
    
    return np.array(celldmsx), np.array(Ex) def read_EtotV(fname):
    """
    Read cell volumes and the corresponding energies from input file *fname*
    (1st col, volumes, 2nd col energies). Units must be :math:`a.u.^3` and 
    :math:`Ryd/cell`
    """
    Vx=[]
    Ex=[]
    
    with open(fname, "r") as lines:
        for line in lines:
            linesplit=line.split()
            V=float(linesplit[0])
            E=float(linesplit[1])
            Vx.append(V)
            Ex.append(E)
    
    return np.array(Vx), np.array(Ex) def read_thermo(fname, ngeo=1):
    """
    Read vibrational thermodynamic functions (Evib, Fvib, Svib, Cvib) as a 
    function of temperature from the input file *fname*. *ngeo* is the number
    of input files to read, corresponding for example to different geometries
    in a quasi-harmonic calculation.
    If *ngeo>1* reads from the files *fname1*, *fname2*, etc. up to *ngeo*  
    Input file(s) have the following format:

    +------+-----------------+-----------------+-----------------+-----------------+
    | T    | :math:`E_{vib}` | :math:`F_{vib}` | :math:`S_{vib}` | :math:`C_{vib}` | 
    +======+=================+=================+=================+=================+
    | 1    | ...             | ...             | ...             | ...             |
    +------+-----------------+-----------------+-----------------+-----------------+  
    
    Lines starting with "#" are not read (comments).    

    Returning values are :math:`nT*ngeo` numpy matrices (T,Evib,Fvib,Svib,Cvib) containing the 
    temperatures and the above mentioned thermodynamic functions as for example:
    Fvib[T,geo] -> Fvib at the temperature *T* for the geometry *geo*
    
    Units must be *K* for temperature, *Ryd/cell* for energies, *Ryd/cell/K* for
    entropy and heat capacity. 
    """
    T=[]
    Evib=[]
    Fvib=[]
    Svib=[]
    Cvib=[]
    for i in range(0,ngeo):
        gT=[]
        gEvib=[]
        gFvib=[]
        gSvib=[]
        gCvib=[]
        if (ngeo>1):
            fname2=fname+str(i+1)
        else:
            fname2=fname
        with open(fname2, "r") as lines:
            for line in lines:
                linesplit=line.split()
                if line!="" and linesplit[0]!="#":
                    gT.append( float(linesplit[0]) )
                    gEvib.append( float(linesplit[1]) )
                    gFvib.append( float(linesplit[2]) )
                    gSvib.append( float(linesplit[3]) )
                    gCvib.append( float(linesplit[4]) )
        T.append(gT)
        Evib.append(gEvib)
        Fvib.append(gFvib)
        Svib.append(gSvib)
        Cvib.append(gCvib)
    
    return np.array(T), np.array(Evib), np.array(Fvib), np.array(Svib), np.array(Cvib) def read_dos_geo(fin,ngeo):
    """
    Read the phonon density of states and energies as in :py:func:`read_dos` from *ngeo* input files
    *fin1*, *fin2*, etc. and store it in two numpy matrices which are returned. 
    """
    
    # read the first file to determine the number of DOS points
    E, dos = read_dos(fin+"1") 
     
    gE = np.zeros((len(E),ngeo))
    gdos = np.zeros((len(E),ngeo))   
    gE[:,0] = E
    gdos[:,0] = dos
    
    # then read the rest of files 
    for i in range(1,ngeo):
        EE, ddos = read_dos(fin+str(i+1)) 
        gE[:,i] = EE
        gdos[:,i] = ddos
    
    return gE, gdos def print_polynomial(a,ibrav=4):
    """
    This function prints the fitted polynomial, either quartic or quadratic
    """
    if ibrav in (1,2,3):                # cubic systems (a,a,a)
        print ("\nFitted polynomial is: \n")
        if (len(a)==3):
            print ("p(x1) = "+str(a[0])+" + "+str(a[1])+" * x1 + "+str(a[2])+" * x1^2\n")
        elif (len(a)==5):
            print ("p(x1) = "+str(a[0])+" + "+str(a[1])+" * x1 + "+str(a[2])+" * x1^2 + "
            +str(a[3])+" * x1^3 + "+str(a[4])+" * x1^4 +\n")
        else:
            print ("Polynomial order not implemented")
    elif ibrav in (4,6,7):              # hexagonal or tetragonal systems (a,a,c)
        print ("\nFitted polynomial is: \n")
        if (len(a)==6):
            print ("p(x1,x2) = "+str(a[0])+" + "+str(a[1])+" * x1 + "+str(a[2])+" * x1^2 + "
            +str(a[3])+" * x2 + "+str(a[4])+" * x2^2 + "+str(a[5])+" *x1*x2\n")
        elif (len(a)==15):
            print ("p(x1,x2) = "+str(a[0])+" + "+str(a[1])+" * x1 + "+str(a[2])+" * x1^2 + "
            +str(a[3])+" * x1^3 + "+str(a[4])+" * x1^4 +\n" 
            +str(a[5])+" *x2 + "+str(a[6])+" *x2^2 + "+str(a[7])+" *x2^3 + "+str(a[8])+" *x2^4 +\n" 
            +str(a[9])+" *x1*x2 + "+str(a[10])+" *x1*x2^2 + "+str(a[11])+" *x1*x2^3 + "+str(a[12])
            +" *x1^2*x2 + "+str(a[13])+" *x1^2*x2^2 +\n"+str(a[14])+" *x1^3*x2\n")
        else:
            print ("Polynomial order not implemented")
    elif ibrav in (8,9,10,11):          # orthorombic systems (a,b,c)
        if (len(a)==10):
            print ("p(x1,x2,x3) = "+str(a[0])+" + "+str(a[1])+" * x1 + "+str(a[2])+" * x1^2 + "
            +str(a[3])+" * x2 + "+str(a[4])+" * x2^2 + "+str(a[5])+" *x1*x2\n")
        elif (len(a)==35):
            print ("p(x1,x2,x3) = "+str(a[0])+" + "+str(a[1])+" * x1 + "+str(a[2])+" * x1^2 + "
            +str(a[3])+" * x1^3 + "+str(a[4])+" * x1^4 +\n" 
            +str(a[5])+" *x2 + "+str(a[6])+" *x2^2 + "+str(a[7])+" *x2^3 + "+str(a[8])+" *x2^4 +\n" 
            +str(a[9])+" *x1*x2 + "+str(a[10])+" *x1*x2^2 + "+str(a[11])+" *x1*x2^3 + "+str(a[12])
            +" *x1^2*x2 + "+str(a[13])+" *x1^2*x2^2 +\n"+str(a[14])+" *x1^3*x2\n"
            +str(a[15])+" *x3"+str(a[16])+" *x3^2"+str(a[17])+" *x3^3"+str(a[18])+" *x3^4\n"
            +str(a[19])+" *x1*x3"+str(a[20])+" *x1*x3^2"+str(a[21])+" *x1*x3^3"+str(a[22])+" *x1^2*x3"
            +str(a[23])+" *x1^2*x3^2"+str(a[24])+" *x1^3*x3"+str(a[25])+" *x2*x3"+str(a[26])+" *x2*x3^2\n"
            +str(a[27])+" *x2*x3^3"+str(a[28])+" *x2^2*x3"+str(a[29])+" *x2^2*x3^2"+str(a[30])+" *x2^3*x3\n"
            +str(a[31])+" *x1*x2*x3"+str(a[32])+" *x1^2*x2*x3"+str(a[33])+" *x1*x2^2x3"+str(a[34])+" *x1*x2*x3^2\n"
            )
        else:
            print ("Polynomial order not implemented")
    else:
        print ("ibrav not implememnted yet") def fit_anis(celldmsx, Ex, ibrav=4, out=False, type="quadratic", ylabel="Etot"):  
    """
    An auxiliary function for handling fitting in the anisotropic case
    """
    if out:
        print (type+" fit")
        if type=="quadratic":
            a, chi = fit_quadratic(celldmsx, Ex, ibrav, out, ylabel)
        elif type=="quartic":
            a, chi = fit_quartic(celldmsx, Ex, ibrav, out, ylabel)
        else:
            print ("Fitting type not implemented")
            return None, None

        if chi!=None:
            print_polynomial(a)
            print ("Chi squared: ",chi,"\n")
            return a, chi 

        return a, None
    
    else:
        if type=="quadratic":
            a, chi = fit_quadratic(celldmsx, Ex, ibrav, False, ylabel)
        elif type=="quartic":
            a, chi = fit_quartic(celldmsx, Ex, ibrav, False, ylabel)
        else:
            return None, None

        if chi!=None:
            return a, chi 

        return a, None def expand_quadratic_to_quartic(a):
    """
    This function gets a vector of coefficients from a quadratic fit and turns it
    into a vector of coeffients as from a quartic fit (extra coeffients are set to zero).
    """
    b=[a[0], a[1], a[2], 0.0, 0.0, 
    a[3], a[4], 0.0, 0.0,
    a[5], 0.0, 0.0, 0.0, 0.0,
    0.0]
    return b def find_min(a,ibrav,type,guess=None,method="BFGS",minoptions={}):
    """
    An auxiliary function for handling the minimum search.    
    """
    if type=="quadratic":
        xmin=find_min_quadratic(a,ibrav,guess,method,minoptions)
        fmin=fquadratic(xmin,a)
    elif type=="quartic":
        xmin=find_min_quartic(a,ibrav,guess,method,minoptions)
        fmin=fquartic(xmin,a)
    else:
        return None
    print ("Minimun "+type+": ",xmin,"\tEnergy at the minimum: "+"{:.20e}".format(fmin)+"\n") 

    return xmin, fmin def fit_Murn(V,E,guess=[0.0,0.0,900/RY_KBAR,1.15],lm_pars={}):
    """
    This is the function for fitting with the Murnaghan EOS as a function of volume only.

    The input variable *V* is an 1D array of volumes, *E* are the corresponding 
    energies (or other analogous quantity to be fitted with the Murnaghan EOS.
    *a*
    
    Note: volumes must be in a.u.^3 and energies in Rydberg.
    
    """
    # reasonable initial guesses for EOS parameters
    if guess[0]==0.0:
        guess[0] = E[len(E) / 2]
    if guess[1]==0.0:
        guess[1] = V[len(V) / 2]

    a, pcov = curve_fit(E_MurnV, V, E, p0=guess, **lm_pars)
    
    chi = 0
    for i in range(0,len(V)):
        chi += (E[i]-E_Murn(V[i],a))**2
    
    return a, pcov, chi def print_eos_data(x,y,a,chi,ylabel="Etot"):   
    """
    Print the data and the fitted results using the Murnaghan EOS. It can be used for
    different fitted quantities using the proper ylabel. ylabel can be "Etot", 
    "Fvib", etc.
    """
    print ("# Murnaghan EOS \t\t chi squared= {:.10e}".format(chi))
    print ("# "+ylabel+"min= {:.10e} Ry".format(a[0])+"\t Vmin= {:.10e} a.u.^3".format(a[1])+"\t B0= {:.10e} kbar".format(a[2]*RY_KBAR)
    +"\t dB0/dV= {:.10e}".format(a[3]))
    print (80*"#")
    print ("# V *a.u.^3)","\t\t",ylabel," (Ry)\t\t",ylabel+"fit"," (Ry)\t\t",ylabel+"-"+ylabel+"fit (Ry)\tP (kbar)")
    for i in range(0,len(y)):
        print ("{:.10e}".format(x[i]),"\t", "{:.10e}".format(y[i])+
        "\t {:.10e}".format(E_Murn(x[i],a))+
        "\t {:.10e}".format(y[i]-E_Murn(x[i],a))+
        "\t {:.10e}".format(P_Murn(x[i],a)*RY_KBAR)) def calculate_fitted_points(V,a):
    """
    Calculates a denser mesh of E(V) points (1000) for plotting.
    """
    Vstep = (V[len(V)-1]-V[0])/1000
    Vdense = np.zeros(1000)
    Edensefitted = np.zeros(1000)
    for i in range(0,1000):
        Vdense[i] = V[0] + Vstep*i 
        Edensefitted[i] = E_Murn(Vdense[i],a)
        
    return Vdense, Edensefitted def compute_beta(TT, minT):
    """
    This function computes the volumetric thermal expansion as a numerical
    derivative of the volume as a function of temperature V(T) given in the
    input array *minT*. This array can obtained
    from the free energy minimization which should be done before.
    """
    grad=np.gradient(np.array(minT))  # numerical derivatives with numpy
    betaT = np.array(grad)  # grad contains the derivatives with respect to T
                                # also convert to np.array format    
    return betaT/minT def compute_beta_splines(TT, minT, splinesoptions={}):
    """
    This function computes the volumetric thermal expansion as a numerical
    derivative of the volume as a function of temperature V(T) given in the
    input array *minT*. This array can obtained
    from the free energy minimization which should be done before.
    This function uses splines for smoothing the derivative.
    """

    betaT = np.zeros(len(TT))

    x = np.array(TT)
    y0 = np.array(minT)

    if (splinesoptions == {}):
        tck0 = interpolate.splrep(x, y0)
    else:
        tck0 = interpolate.splrep(x, y0, k=splinesoptions['k0'], s=splinesoptions['s0'])

    ynew0 = interpolate.splev(x, tck0, der=0)
    betaT = interpolate.splev(x, tck0, der=1)

    return betaT/minT def compute_Cv(T,Vmin,V,Cvib):
    """
    This function computes the isocoric heat capacity as a function of temperature.
    From *Cvib*, which is a matrix with *Cvib(T,V)* as from the harmonic calculations
    determines the *Cv* at each temperature by linear interpolation between the values
    at the two volumes closest to Vmin(T). Vmin(T) is from the minimization of F(V,T)
    and *V* is the array of volumes used for it.
    Returns *Cv(T)*.
    
    Work in progress... for now it uses all volumes in the interpolation.
    
    """
    Cv = np.zeros(len(T))
    for iT in range(0,len(T)):
        Cv_interpolated = np.interp(Vmin[iT], V, Cvib[iT,:])
        Cv[iT] = Cv_interpolated
        
    return Cv def compute_Cp(T,Cv,V,B0,beta):
    """
    This function computes the isobaric heat capacity from the equation:
    
    :math:`Cp - Cv = T V beta^2 B0`

    where *Cp,Cv* are the isobaric and isocoric heat capacities respectively,
    *T* is the temperature, *V* the unit cell volume, *beta* the volumetric
    thermal expansion and *B0* the isothermal bulk modulus.
    
    """
    Cp = Cv + T * V * beta * beta * B0
    return Cp def compute_alpha(minT,ibrav):
    """
    This function calculates the thermal expansion alphaT at different temperatures
    from the input minT matrix by computing the numerical derivatives with numpy.
    The input matrix minT has shape nT*6, where the first index is the temperature 
    and the second the lattice parameter. For example, minT[i,0] and minT[i,2] are
    the lattice parameters a and c at the temperature i.
 
    More ibrav types must be implemented
    """
    
    grad=np.gradient(np.array(minT))  # numerical derivatives with numpy
    alphaT = np.array(grad[0])  # grad[0] contains the derivatives with respect to T, which is the first axis in minT
                                # also convert to np.array format
    
    # now normalize the alpha properly. It must be different for different ibrav
    # to avoid a divide by 0 error (minT is zero for lattice parameters not defined
    # in the system)
    if ibrav==4:
        alphaT[:,0] = alphaT[:,0]/minT[:,0]
        alphaT[:,2] = alphaT[:,2]/minT[:,2]
        
    return alphaT def compute_alpha_splines(TT,minT,ibrav,splinesoptions):
    """
    This function calculates the thermal expansions alphaT at different temperatures
    as the previous function but using spline interpolation as implemented in
    scipy.interpolate.

    """
    alphaT = np.zeros(len(TT)*6)
    alphaT.shape = (len(TT),6)
    
    x = np.array(TT)
    y0 = np.array(minT[:,0])
    y1 = np.array(minT[:,1])
    y2 = np.array(minT[:,2])
        
    if (splinesoptions=={}):
        tck0 = interpolate.splrep(x, y0)
        tck1 = interpolate.splrep(x, y1)
        tck2 = interpolate.splrep(x, y2)  
    else:
        tck0 = interpolate.splrep(x, y0, k=splinesoptions['k0'], s=splinesoptions['s0'])
        tck1 = interpolate.splrep(x, y1, k=splinesoptions['k1'], s=splinesoptions['s1'])
        tck2 = interpolate.splrep(x, y2, k=splinesoptions['k2'], s=splinesoptions['s2'])        
        
    ynew0 = interpolate.splev(x, tck0, der=0)
    alphaT[:,0] = interpolate.splev(x, tck0, der=1)
    ynew1 = interpolate.splev(x, tck1, der=0)
    alphaT[:,1] = interpolate.splev(x, tck1, der=1)
    ynew2 = interpolate.splev(x, tck2, der=0)
    alphaT[:,2] = interpolate.splev(x, tck2, der=1)
    
    # now normalize the alphaTs properly. It must be different for different ibrav
    # to avoid a divide by 0 error (minT is zero for lattice parameters not defined
    # in the system)
    if ibrav==4:
        alphaT[:,0] = alphaT[:,0]/minT[:,0]
        alphaT[:,2] = alphaT[:,2]/minT[:,2]
        
    return alphaT def write_celldmsT(fname,T,x,ibrav=4):
    fout=open(fname,"w")
    if ibrav==4:    # Hexagonal systems    
        fout.write("# T (K)\t\ta\t\tc\n") 
        for i in range(0,len(T)): 
            fout.write("{:.10e}".format(T[i])+"\t"+"{:.10e}".format(x[i][0])+"\t"+
            "{:.10e}".format(x[i][2])+"\t"+"\n") 
    fout.close() def write_alphaT(fname,T,alphaT,ibrav=4):
    fout=open(fname,"w")
    if ibrav==4:    # Hexagonal systems    
        fout.write("# T (K) \t\talpha_xx\t\talpha_zz\n") 
        for i in range(0,len(T)): 
            fout.write("{:.10e}".format(T[i])+"\t"
            +"\t"+"{:.10e}".format(alphaT[i,0])+"\t"+"{:.10e}".format(alphaT[i,2])+"\n") 
    fout.close() def write_xy(fname,x,y,labelx,labely):
    """
    This function writes a quantity y versus quantity x into the file fname.
    y and x are arrays and should have the same lenght. labelx and labely
    are the axis labels (possibly with units), written in the header of the file
    (first line).
    """
    fout=open(fname,"w")
    header = "# "+labelx+"\t"+labely+"\n"
    fout.write(header)     
    for i in range(0,len(x)): 
        line = "{:.10e}".format(x[i])+"\t\t"+"{:.10e}".format(y[i])+"\n"
        fout.write(line)
    fout.close() def gen_TT(Tstart=1,Tend=1000,Tstep=1):
    """
    A simple function to generate a numpy array of temperatures, starting from
    *Tstart* and ending to *Tend* (or the closest *T<Tend* accorinding to the *Tstep* )
    with step *Tstep* .
    """
    Tsteps = int((Tend-Tstart)//Tstep)
    TT = np.zeros(Tsteps)
    for i in range(0,Tsteps):
        TT[i] = Tstart + i*Tstep
        
    return TT def compute_thermo_geo(fin,fout=None,ngeo=1,TT=np.array([1])):
    """
    This function reads the input dos file(s) from *fin+i*, with *i* a number from
    1 to *ngeo* + 1 and computes vibrational energy, Helmholtz energy, entropy and
    heat capacity in the harmonic approximation. Then writes the output on file(s)
    if fout!=None.
    Output file(s) have the following format:

    +------+-----------------+-----------------+-----------------+-----------------+
    | T    | :math:`E_{vib}` | :math:`F_{vib}` | :math:`S_{vib}` | :math:`C_{vib}` | 
    +======+=================+=================+=================+=================+
    | 1    | ...             | ...             | ...             | ...             |
    +------+-----------------+-----------------+-----------------+-----------------+   

    and are names *fout* +1, *fout* +2,... for each geometry.

    Returning values are (len(TT),ngeo) numpy matrices (T,gEvib,gFvib,gSvib,gCvib,gZPE,gmodes) 
    containing the 
    temperatures and the above mentioned thermodynamic functions as for example:
    Fvib[T,geo] -> Fvib at the temperature "T" for the geometry "geo"
    """

    gEvib=np.zeros((len(TT),ngeo))
    gFvib=np.zeros((len(TT),ngeo))
    gSvib=np.zeros((len(TT),ngeo))
    gCvib=np.zeros((len(TT),ngeo))
    gZPE=np.zeros((ngeo))
    gmodes=np.zeros((ngeo))
    for i in range(0,ngeo):
        E, dos = read_dos(fin+str(i+1))        
        T, Evib, Svib, Cvib, Fvib, ZPE, modes = compute_thermo(E/RY_TO_CMM1,dos*RY_TO_CMM1,TT)
        if (fout!=None):
            write_thermo(fout+str(i+1),T, Evib, Fvib, Svib, Cvib, ZPE, modes)
        
        gEvib[:,i]=Evib
        gFvib[:,i]=Fvib
        gSvib[:,i]=Svib
        gCvib[:,i]=Cvib
        gZPE[i]=ZPE
        gmodes[i]=modes
    
    return TT, gEvib, gFvib, gSvib, gCvib, gZPE, gmodes def rearrange_thermo(T,Evib,Fvib,Svib,Cvib,ngeo=1):
    """
    This function just rearranges the order of the elements in the input matrices
    The first index of the returning matrices *X* now gives all geometries at a given
    *T*, i.e. *X[0]* is the vector of the property *X* a *T=T[0,0]* . *X[0,0]* for the first 
    geometry, *X[0,1]* the second geometry and so on.
    """
    Evib2 = np.zeros((len(T[0]),ngeo))
    Fvib2 = np.zeros((len(T[0]),ngeo))
    Svib2 = np.zeros((len(T[0]),ngeo))
    Cvib2 =np.zeros((len(T[0]),ngeo)) 

    for i in range(0,len(T[0])):
        for j in range(0,ngeo):
            Evib2 [i,j] = Evib[j][i]
            Fvib2 [i,j] = Fvib[j][i]
            Svib2 [i,j] = Svib[j][i]
            Cvib2 [i,j] = Cvib[j][i]
    
    return len(T[0]), T[0], Evib2, Fvib2, Svib2, Cvib2 def simple_plot_xy(x,y,xlabel="",ylabel=""):
    """
    This function generates a simple xy plot with matplotlib.
    """
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1)  # create an axes object in the figure
    ax.plot(x, y, 'r')
    ax.set_xlabel(xlabel)
    ax.set_ylabel(ylabel)
    plt.show()
    
    return fig def multiple_plot_xy(x,y,xlabel="",ylabel="",labels=""):
    """
    This function generates a simple xy plot with matplotlib overlapping several
    lines as in the matrix y. y second index refers to a line in the plot, the first 
    index is for the array to be plotted.
    """
    
    if (len(y[0,:])>7):
        print ("Too many data on y axis!")
        return
    
    colors = ['k','r','b','g','c','m','y']
    
    fig = plt.figure()
    ax = fig.add_subplot(1, 1, 1)  # create an axes object in the figure
    if (labels==""):     
        try:    # try if there are multiple data on x axis
            for i in range(0,len(y[0,:])):
                ax.plot(x[:,i], y[:,i], colors[i])
        except: # if not use a single x axis
            for i in range(0,len(y[0,:])):
                ax.plot(x, y[:,i], colors[i])
    else:
        try:    # try if there are multiple data on x axis
            for i in range(0,len(y[0,:])):
                ax.plot(x[:,i], y[:,i], colors[i],label=labels[i])
        except: # if not use a single x axis
            for i in range(0,len(y[0,:])):
                ax.plot(x, y[:,i], colors[i],label=labels[i])       
        ax.legend() 
        
    ax.set_xlabel(xlabel)
    ax.set_ylabel(ylabel)

    plt.show()
    
    return fig
import sys
import numpy as np
from .constants import RY_KBAR, RYDBERG_SI, NA
from .read import read_Etot, read_EtotV, read_thermo, read_dos_geo
from .fitutils import print_polynomial, fit_anis, expand_quadratic_to_quartic
from .minutils import find_min
from .eos import fit_Murn, print_eos_data, calculate_fitted_points, compute_beta, compute_beta_splines, compute_Cv, compute_Cp
from .properties_anis import compute_alpha, compute_alpha_splines
from .write import write_celldmsT, write_alphaT, write_xy
from .thermo import gen_TT, compute_thermo_geo, rearrange_thermo
from .plotutils import simple_plot_xy, multiple_plot_xy
    *k1*, *s1*, *k2*, *s2* which are passed to :py:func:`scipy.interpolate.splrep`, 
    one couple for each set of thermal expansions (alpha_xx, alpha_yy, alpha_zz).
    *k* is the order of the spline (default=3), s a smoothing condition (default=None).
    If *splinesoptions=={}* use the default options of :py:func:`scipy.interpolate.splrep`
    Note: use this option with care
       
    """
        
    print ("Minimization method: ",method,"\t options: ",minoptions)
    
    # Read the Etot at 0 K
    celldmsx, Ex = read_Etot(fEtot)
    
    ngeo = len(Ex)   # total number of geometries

    nT, T, Evib, Fvib, Svib, Cvib = thermodata
    
    # Fit and find the minimun at 0 K
    a0temp, chi0 = fit_anis(celldmsx, Ex, ibrav, out=True, type=typeEtot, ylabel="Etot")
    if chi0!=None:
        min0, fmin0 = find_min(a0temp, ibrav, type=typeEtot, guess=defaultguess, method=method,minoptions=minoptions)
     
    # if a mixed type fitting has been used, expand the vector of polynomial coefficients
    if (typeEtot=="quadratic" and typeFvib=="quartic"):
        a0=expand_quadratic_to_quartic(a0temp)
    else:
        a0=a0temp
        
    # Fit and find the minimun at T>0 K for each T 
    minT = np.zeros((nT,6))
    fminT= np.zeros(nT)
    aT = np.zeros((nT,len(a0)))
    chiT= np.zeros(nT)
    minguess = min0
    for i in range(0,nT): 
        print ("####################################################################")
        print ("T=",T[i])
  
        # Fit Fvib with a quadratic or quartic polynomial, then add the fitted 
        # quadratic or quartic Etot at 0 K and find the minimun
        aTtemp, chiT[i] = fit_anis(celldmsx,Fvib[i],ibrav,type=typeFvib, ylabel="Fvib")
        if (chiT[i]>0):
            print (typeEtot+" Etot + "+typeFvib+" Fvib polynomials")
            if (typeEtot=="quartic" and typeFvib=="quadratic"):
                aT[i,:]=expand_quadratic_to_quartic(aTtemp)
                typemix="quartic"       # find the minimum as a quadratic
            else:
                aT[i,:]=aTtemp       
                typemix=typeFvib         # typeFvib and typeFtot are the same
            print_polynomial(a0+aT[i,:])

            if (T[i]>50):   # Tests have shown that for temperatures > 50~K  you get better results 
                            # using as initial guess the minimum at the previous temperature
                minguess=minT[i-1,:]
            minT[i,:], fminT[i] = find_min(a0+aT[i,:],ibrav,type=typemix,guess=minguess,method=method,minoptions=minoptions)            

    # Calculate the thermal expansion(s)
    if (splinesoptions==None):
        alphaT = compute_alpha(minT,ibrav)
    else:
        alphaT = compute_alpha_splines(T,minT,ibrav,splinesoptions)
    
    return T, fminT, minT, alphaT, a0, chi0, aT, chiT


################################################################################

def fitFvibV(fin,thermodata,splinesoptions={},lm_pars={},verbosity="low"):
    """
    This function computes quasi-harmonic quantities from the 
    :math:`E_{tot}(V)+F_{vib}(V,T)` as a function of temperature with Murnaghan's
    EOS. :math:`E_{tot}(V)` is read from the *fin* file. :math:`F_{vib}(V,T)`
    are given in *thermodata* which is a list containing the number of temperatures
    ( *nT* ) for which the calculations are done and the numpy matrices for 
    temperatures, vibrational energy, Helmholtz energy, entropy and
    heat capacity. All these quantities are for each volume as in *fin* file.
        
    The function fits :math:`E_{tot}(V)+F_{vib}(V,T)` with a Murnaghan's EOS
    at each temperature in *thermodata* and then stores the fitted coefficients.
    It also computes the volume thermal expansion as a numerical derivative of
    the minimum volume as a function of temperature (:py:func:`compute_beta`), the
    constant volume heat capacity at the minimum volume at each T
    (:py:func:`compute_Cv`) and the constant pression heat capacity (:py:func:`compute_Cp`).
    
    It returns the numpy 1D arrays containing the temperatures (as in input), the
    minimun energy,lm_pars minimun volume, bulk modulus, volume thermal expansion, constant
    volume and constant pressure heat capacities, one matrix with all fitted 
    coefficients at each T and finally an array with the :math:`\chi^2` at each T.
    
    .. Warning::
       The quantities in *thermodata* are usually obtained from :py:func:`compute_thermo_geo`
       or from :py:func:`read_thermo` and :py:func:`rearrange_thermo`. It is
       important that the order in the total energy file *fin* and the order of
       the thermodynamic data in *thermodata* is the same!  See also *example5* and 
       the tutorial. 
       
    """
    
    V, E = read_EtotV(fin)
    
    a0, cov, chi = fit_Murn(V,E)
   
    print_eos_data(V,E,a0,chi,"E")

    ngeo = len(E)   # total number of geometries
    
    nT, T, Evib, Fvib, Svib, Cvib = thermodata  # import list of thermodynamic data
   
    # Fit and find the minimun at T>0 K for each T                    
    aT = np.zeros((nT,4))    # Murnaghan EOS always has 4 parameters, store one for eact T
    chiT = np.zeros((nT))
    for i in range(0,nT): 
        print ("####################################################################")
        print ("T=",T[i])
  
        a, cov, chi = fit_Murn(V,E+Fvib[i],a0,lm_pars)
        if (verbosity=="high"):
            print_eos_data(V,E+Fvib[i],a,chi,"E+Fvib")  # print full detail at each T
             
        print ("Fmin = {:.10e}".format(a[0])+"\tVmin = {:.10e}".format(a[1])

        +"\tB0 = {:.10e}".format(a[2]*RY_KBAR)+"\t chi^2= {:.10e}".format(chi))