# Path: pyqha/constants.py
# RY_TO_CMM1 = 1.0E+10 * RY_TO_THZ / C_SI
# 
# K_BOLTZMANN_RY   = K_BOLTZMANN_SI / RYDBERG_SI
# 
# Path: pyqha/read.py
# def read_dos(filename):
#     """
#     Read the phonon density of states (y axis) and the corresponding energies (x axis)
#     from the input file *filename* (1st col energies, 2nd col DOS) and store it
#     in two numpy arrays which are returned. 
# 
#     """
#     E = []
#     dens = []
#     print ("Reading phonon dos file "+filename+"...")
#     with open(filename, "r") as lines:
#         for line in lines:
#             linesplit=line.split()
#             if (float(linesplit[0])<1.0E-10):
#                 E.append(1.0E-10)
#             else:
#                 E.append(float(linesplit[0]))
#             dens.append(float(linesplit[1]))
# 
#     return np.array(E), np.array(dens)
# 
# 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
# 
# Path: pyqha/write.py
# def write_thermo(fname,T, Evib, Fvib, Svib, Cvib,ZPE,modes):
#     fout=open(fname,"w")
#     fout.write('# total modes from dos = {:.10e}\n'.format(modes))
#     fout.write('# ZPE = {:.10e} Ry/cell\n'.format(ZPE))
#     fout.write("# Multiply by 13.6058 to have energies in eV/cell etc..\n\
# # Multiply by 13.6058 x 23060.35 = 313 754.5 to have energies in cal/(N mol).\n\
# # Multiply by 13.6058 x 96526.0 = 1 313 313 to have energies in J/(N mol).\n\
# # N is the number of formula units per cell.\n#\n")
#     
#     fout.write("# T (K) \tEvib (Ry/cell)\tFvib (Ry/cell)\tSvib (Ry/cell/K)\tCvib (Ry/cell/K)\n") 
#     for i in range(0,len(Evib)): 
#         fout.write('{:.10e}\t{:.10e}\t{:.10e}\t{:.10e}\t{:.10e}\n'.format(T[i],Evib[i],Fvib[i],Svib[i],Cvib[i]))
#     fout.close()


from .constants import RY_TO_CMM1, K_BOLTZMANN_RY
from math import tanh, sinh, log, exp
from .read import read_dos, read_dos_geo
from .write import write_thermo
import numpy as np
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))        

# Path: pyqha/constants.py
# RY_TO_CMM1 = 1.0E+10 * RY_TO_THZ / C_SI
# 
# K_BOLTZMANN_RY   = K_BOLTZMANN_SI / RYDBERG_SI
# 
# Path: pyqha/read.py
# def read_dos(filename):
#     """
#     Read the phonon density of states (y axis) and the corresponding energies (x axis)
#     from the input file *filename* (1st col energies, 2nd col DOS) and store it
#     in two numpy arrays which are returned. 
# 
#     """
#     E = []
#     dens = []
#     print ("Reading phonon dos file "+filename+"...")
#     with open(filename, "r") as lines:
#         for line in lines:
#             linesplit=line.split()
#             if (float(linesplit[0])<1.0E-10):
#                 E.append(1.0E-10)
#             else:
#                 E.append(float(linesplit[0]))
#             dens.append(float(linesplit[1]))
# 
#     return np.array(E), np.array(dens)
# 
# 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
# 
# Path: pyqha/write.py
# def write_thermo(fname,T, Evib, Fvib, Svib, Cvib,ZPE,modes):
#     fout=open(fname,"w")
#     fout.write('# total modes from dos = {:.10e}\n'.format(modes))
#     fout.write('# ZPE = {:.10e} Ry/cell\n'.format(ZPE))
#     fout.write("# Multiply by 13.6058 to have energies in eV/cell etc..\n\
# # Multiply by 13.6058 x 23060.35 = 313 754.5 to have energies in cal/(N mol).\n\
# # Multiply by 13.6058 x 96526.0 = 1 313 313 to have energies in J/(N mol).\n\
# # N is the number of formula units per cell.\n#\n")
#     
#     fout.write("# T (K) \tEvib (Ry/cell)\tFvib (Ry/cell)\tSvib (Ry/cell/K)\tCvib (Ry/cell/K)\n") 
#     for i in range(0,len(Evib)): 
#         fout.write('{:.10e}\t{:.10e}\t{:.10e}\t{:.10e}\t{:.10e}\n'.format(T[i],Evib[i],Fvib[i],Svib[i],Cvib[i]))
#     fout.close()


# Path: pyqha/thermo.py
from .constants import RY_TO_CMM1, K_BOLTZMANN_RY
from math import tanh, sinh, log, exp
from .read import read_dos, read_dos_geo
from .write import write_thermo
import numpy as np

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)