import CDF
import numpy as np
import matplotlib.pyplot as plt
def PlotDFs(DFs,hardcopy=False,title='',linetype='',nolegendintra=False,figsize=(20,14),dpi=80):
    '''Plot all objects of class CDF (distribution functions: RDF,potential, reference RDF, corrected potential) presented in list DFs on same plot,
    Hardcopy - create eps-copy of the plot,
    linetype - define type of linse using same syntax as matplotlib.pyplot
    nolegendintra - turn off legend on plots with intramolecular DFs. Usually it is usefull then legend overlaps with lines on plot. 
    figsize, dpi - sizeand resolution of the resulting plot in inches and dpi '''
    fig=plt.figure(figsize=figsize,dpi=dpi)
    if isinstance(DFs,type([])):
        plt.hold(False)        
        for i in DFs:
            if linetype=='':
                i.Plot(hidden=True,nolegendintra=nolegendintra)
            else:
                i.Plot(hidden=True,linetype=linetype,nolegendintra=nolegendintra)
            plt.hold(True)            
        plt.title(title)
        plt.xlabel('r, A')
        plt.legend(loc=0)
        if hardcopy:
             plt.savefig(title.replace(' ','.')+'.eps')
        plt.show()
    if isinstance(DFs,CDF.C_DF):
        DFs.Plot(hardcopy=hardcopy)
def AnalyzeIMCOuput(filename,DF='RDF',mtype='dmpc',iters=()):
    '''Analyze output file <filename> produced by MAGIC software and plot resulting DFs of interest
    DF - defines function of interest: RDF-radial (and angle) distribution functions obtained in inverse MCcalculation ,RDFref-reference RDFs(ADFs), Pot - Potential used in a current iteration, PotNew - potential generated in current iteration.
    mtype - type of the system. It defines list of pairs of atomic names to search in output file. All available values of mtype can be found in help of GetDFsFromFile_magic
    iters() - tuple of iterations to extract from the file. If nothing mentioned all iterations will be extracted.  Example: iters=(1,2,3) or iters=(1,) Mind the comma after single iteration number!
    '''
    if DF.upper()=='RDF':
        DFs=GetDFsFromFile_magic(filename,iters=iters,DF='RDF',mtype=mtype,basename=filename)
        DFsRef=GetDFsFromFile_magic(filename,iters=iters,DF='RDFref',mtype=mtype,basename=filename)
        for i in xrange(0,len(DFs[0])):
            lst=[DFs[j][i] for j in xrange(0,len(DFs))]
            lst.insert(0,DFsRef[0][i])
            PlotDFs(lst,hardcopy=True)
    
def DumpDFs(DFs,filename):
    '''Dump a list of distribution functions DFs to a file <filename> using module "pickle".'''
    import pickle
    f=open(filename,'ab')
    pickle.dump(DFs,f)
    f.close()
def LoadDFs(filename):
    '''Load a list of distribution functions DFs from a file <filename> using module "pickle".'''
    import pickle
    f=open(filename,'rb')
    DFs=pickle.load(f)
    f.close()
    return DFs
def GetOptEpsilon(pots,eps_old,r1,eps_min=0,eps_max=0,npoints=100):
    '''
    Calculate optimal value of dielectric permittivity which provides fastest decay of short-range intermolecular potentials.
    Pots - set of potentials
    Eps_old is a dielectric permittivity used for IMC calcuation
    R1 - range defining where potential tail starts.
    eps_min,eps_max - define range of eps search, by default eps_min=0, eps_max=2*eps_old
    npoints - how many points in the eps_mesh to have
    '''
    
    # Generate mesh of epsilon values:
    weight=1e25
    eps_opt=0
    if (eps_min==0) and(eps_max==0):
        eps_max=eps_old*2
    eps=np.linspace(eps_min,eps_max,npoints)
    for e in eps: #every particular value of the eps-mesh
        w=0.0
        for pot in pots: # loop over all potentials 
            if pot.type=='inter':
                w+=pot.GetWeight(e,float(eps_old),r1)
        if w<weight:
            weight=w
            eps_opt=e
#        print 'eps: %5.3f  weight: %12.7f' %(e,w)
    return eps_opt
def PotsEpsCorrection(pots,eps_old,eps_new):
    '''
    Create new short-range potentials with corrected value of dielectric permittivity for all intermolecular potentials in list "pots". Intramolecular potentials are kept untouched
    '''
    import copy
    newpots=copy.deepcopy(pots)
    for pot in newpots:
        if pot.type=='inter':
            U_sr_new=pot.g[:,1]+((pot.qq*1.602e-19*1.602e-19*6.022e23)/(4.0*3.1416*8.85e-12*1.0e-10*1.0e3))*(1.0/eps_old-1.0/eps_new)*(1.0/pot.g[:,0])
            pot.g[:,1]=U_sr_new
            pot.Name=pot.Name+' eps_old:%5.3f eps_new: %5.3f' %(eps_old,eps_new)
    return newpots
def TotalPots(pots,eps):
    '''
    Create a set ot total potentials by including electrostatics to short-range interatomic potentials
    '''
    import copy
    newpots=copy.deepcopy(pots)
    for pot in newpots:
        if pot.type=='inter':
            U_sr_new=pot.g[:,1]+((pot.qq*1.602e-19*1.602e-19*6.022e23)/(4.0*3.1416*8.85e-12*1.0e-10*1.0e3))*(1.0/eps)*(1.0/pot.g[:,0])
            pot.g[:,1]=U_sr_new
            pot.Name=pot.Name+' Total with eps: %5.3f' %(eps)
    return newpots

def GetDFsFromFile_rdf(rdffilename,basename=''):
    '''
    Analyse rdf-file <rdffilename> produced by rdf-2.0 or rdf-1.3 tool and return a list of DFs.
    basename - prefix given to all the DFs readed from the file
    '''
    RDFsList=[]
    rdffile=open(rdffilename,'r')
    lines=rdffile.readlines()
    line=lines.pop(0)
    NAtomTypes=int(line.split()[0])
    NPoints=int(line.split()[1])
    RMin=float(line.split()[2])
    RMax=float(line.split()[3])
    if NAtomTypes>0:
        resol=(RMax-RMin)/NPoints
    AtomTypes=[]
    RDFs=np.zeros((NAtomTypes,NAtomTypes,NPoints,2))
    for i in xrange(0,NAtomTypes):
        line=lines.pop(0).split('\n')[0]
        AtomTypes.append(line.strip())
    while (lines!=[]):
        line=lines.pop(0).split('\n')[0]
        if (line.strip().upper()=='INTRA') or (line.strip().upper()=='ANGLE') or (line.strip().upper()=='END'):
             break
        r=float(line.split()[0])
        g_r=float(line.split()[1])
        i=int(line.split()[2])-1
        j=int(line.split()[3])-1
        ind=int((r-RMin)/resol)
        RDFs[i,j,ind,0]=r
        RDFs[i,j,ind,1]=g_r     
    for i in xrange(0,NAtomTypes):
        for j in xrange(i,NAtomTypes):
            trdf=CDF.C_DF('inter',basename+'RDF '+AtomTypes[i]+'-'+AtomTypes[j],RMin,RMax,NPoints,g=RDFs[i,j,:,:])
            RDFsList.append(trdf)
    if (line.strip().upper()=='INTRA'):
        while (lines!=[]):
            line=lines.pop(0).split('\n')[0]
            if (line.strip().upper()=='ANGLE') or (line.strip().upper()=='END'):
                break
            npairs=int(line.split()[0])
            NPoints=int(line.split()[1])
            RMin=float(line.split()[2])
            RMax=float(line.split()[3])
            resol=(RMax-RMin)/NPoints
            RDF=np.zeros((NPoints,2))
            name=''
            for i in xrange(0,npairs):
                line=lines.pop(0).split('\n')[0]
            for i in xrange(0,NPoints):
                line=lines.pop(0).split('\n')[0]
                r=float(line.split()[0])
                g_r=float(line.split()[1])
                if name=='':
                    for k in xrange(2,len(line.replace('rdf_int:','').split())):
                        name=name+line.replace('rdf_int:','').split()[k].strip()
                ind=int((r-RMin)/resol)
                RDF[ind,0]=r
                RDF[ind,1]=g_r
            trdf=CDF.C_DF('intra',basename+'RDF '+name,RMin,RMax,NPoints,g=RDF)
            RDFsList.append(trdf)
            
    if (line.strip().upper()=='ANGLE'):
        while (lines!=[]):
            line=lines.pop(0).split('\n')[0]
            if (line.strip().upper()=='ANGLE') or (line.strip().upper()=='END'):
                break
            npairs=int(line.split()[0])
            NPoints=int(line.split()[1])
            RMin=float(line.split()[2])
            RMax=float(line.split()[3])
            resol=(RMax-RMin)/NPoints
            RDF=np.zeros((NPoints,2))
            name=''
            for i in xrange(0,npairs):
                line=lines.pop(0).split('\n')[0]
            for i in xrange(0,NPoints):
                line=lines.pop(0).split('\n')[0]
                r=float(line.split()[0])
                g_r=float(line.split()[1])
                if name=='':
                    for k in xrange(2,len(line.split())):
                        name=name+line.split()[k].strip()
                ind=int((r-RMin)/resol)
                RDF[ind,0]=r
                RDF[ind,1]=g_r
            trdf=CDF.C_DF('intra',basename+'ADF '+name,RMin,RMax,NPoints,g=RDF)
            RDFsList.append(trdf)
    return RDFsList
def GetDFsFromFile_dat(rdffilename,basename=''):
    '''
    Analyse dat-file <rdffilename> produced by  rdf-1.3 tool and return a list of DFs.
    basename - prefix given to all the DFs readed from the file
    '''
    import re
    RDFsList=[]
    rdffile=open(rdffilename,'r')
    lines=rdffile.readlines()
    PairNamesList=[]
    for l in lines:
        if re.search('This pair:',l):
            PairNamesList.append(l.split('\n')[0].split('This pair:')[1])

#Loop over intermolecular RDFs
    for PairName in PairNamesList:                                             
        if re.search('Npairs',PairName):
            RDF_r_List=[]
            RDF_gr_List=[]
            sstring='rdf_smooth: *'+PairName.split('Npairs')[0].strip().split('-')[0].strip()+' *- *'+PairName.split('Npairs')[0].strip().split('-')[1].strip()
            name=PairName.split('Npairs')[0].strip().split('-')[0].strip()+'-'+PairName.split('Npairs')[0].strip().split('-')[1].strip()
            for l in lines:
                if re.search(sstring,l):
                    RDF_r_List.append(float(l.split()[0]))
                    RDF_gr_List.append(float(l.split()[1]))
            RDF=np.array([[RDF_r_List[i],RDF_gr_List[i]] for i in xrange(0,len(RDF_r_List))])
            RMin=RDF[:,0].min()
            RMax=RDF[:,0].max()
            Npoints=len(RDF[:,0])
            resol=(RMax-RMin)/(Npoints-1)
            for i in xrange(0,Npoints-1):
                 if (abs(RDF[i+1,0]-RDF[i,0]-resol)>=0.0000001):
                     print 'Error in RDF, points are not unidistant'
                     return
            trdf=CDF.C_DF('inter',basename+'RDF '+name,RMin,RMax,Npoints,g=RDF)
            RDFsList.append(trdf)
        if re.search('Nint',PairName):
            RDF_r_List=[]
            RDF_gr_List=[]
            sstring='rdf_int_smooth: *'+PairName.split('Nint')[0].strip().split('-')[0].strip()+' *- *'+PairName.split('Nint')[0].strip().split('-')[1].strip()
            name=PairName.split('Nint')[0].strip().split('-')[0].strip()+'-'+PairName.split('Nint')[0].strip().split('-')[1].strip()
            for l in lines:
                if re.search(sstring,l):
                    RDF_r_List.append(float(l.split()[0]))
                    RDF_gr_List.append(float(l.split()[1]))
            RDF=np.array([[RDF_r_List[i],RDF_gr_List[i]] for i in xrange(0,len(RDF_r_List))])
            RMin=RDF[:,0].min()
            RMax=RDF[:,0].max()
            Npoints=len(RDF[:,0])
            resol=(RMax-RMin)/(Npoints-1)
            RMin=RMin-resol/2
            RMax=RMax+resol/2
            for i in xrange(0,Npoints-1):
                if (abs(RDF[i+1,0]-RDF[i,0]-resol)>=0.0000001):
                     print 'Error in RDF, points are not unidistant'
                     return
            trdf=CDF.C_DF('intra',basename+'RDF '+name,RMin,RMax,Npoints,g=RDF)
            RDFsList.append(trdf)
    return RDFsList

def GetDFsFromFile_magic(rdffilename,iters=(),DF='RDF',mtype='dmpc',basename=''):
    '''
    Analyse file <rdffilename> produced by MAGIC software and return a list of lists of DFs. DFs[iteration][pair]
    iters() - tuple of iterations to extract from the file. If nothing mentioned all iterations will be extracted.  Example: iters=(1,2,3) or iters=(1,) Mind the comma after single iteration number!
    basename - prefix given to all the DFs readed from the file    
    DF - defines function of interest: RDF-radial (and angle) distribution functions obtained in inverse MCcalculation ,RDFref-reference RDFs(ADFs), Pot - Potential used in a current iteration, PotNew - potential generated in current iteration.
    mtype - type of the system. It defines list of pairs of atomic names to search in output file. 
        the list consists of three sublists: first of them refer to non-bonded pair interactions, the second refers to bonded pair interactions, and the third refers to bending angle (1-3) bond interactions
        if mtype=='dmpc':
            PairNamesList=[['N-N','N-P','N-CH','N-CO','P-P','P-CH','P-CO','CH-CH','CH-CO','CO-CO'],['N-P','P-C1','C2-C3','C1-C2','C1-C5'],['N-C1-P','P-C2-C1','C1-C3-C2','C2-C4-C3']]    
        if mtype=='dmpc_12': No 1-3 interaction
            PairNamesList=[['N-N','N-P','N-CH','N-CO','P-P','P-CH','P-CO','CH-CH','CH-CO','CO-CO'],['N-P','P-C1','C2-C3','C1-C2','C1-C5'],[]]    
    '''
    import re
    if mtype=='dmpc':
        PairNamesList=[['N-N','N-P','N-CH','N-CO','P-P','P-CH','P-CO','CH-CH','CH-CO','CO-CO'],['N-P','P-C1','C2-C3','C1-C2','C1-C5'],['N-C1-P','P-C2-C1','C1-C3-C2','C2-C4-C3']]    
    if mtype=='hxd':
        PairNamesList=[['CH-CH'],['C1-C2'],[]]    
    if mtype=='hxd_A':
        PairNamesList=[['CH-CH'],['C1-C2'],['C1-C3-C2']]    
    if mtype=='dmpc_12':
        PairNamesList=[['N-N','N-P','N-CH','N-CO','P-P','P-CH','P-CO','CH-CH','CH-CO','CO-CO'],['N-P','P-C1','C2-C3','C1-C2','C1-C5'],[]]    
    if mtype=='dmpc_A6':
        PairNamesList=[['N-N','N-P','N-CH','N-CO','P-P','P-CH','P-CO','CH-CH','CH-CO','CO-CO'],['N-P','P-C1','C2-C3','C1-C2','C1-C5'],['N-C1-P','P-C2-C1','C1-C3-C2','C2-C4-C3','C2-C5-C1','C6-C1-C5']]    
    if mtype=='dmpc_A5':
        PairNamesList=[['N-N','N-P','N-CH','N-CO','P-P','P-CH','P-CO','CH-CH','CH-CO','CO-CO'],['N-P','P-C1','C2-C3','C1-C2','C1-C5'],['N-C1-P','P-C2-C1','C1-C3-C2','C2-C4-C3','C2-C5-C1']]    

    if basename!='':
        basename+=' '
    if DF.upper() not in ('RDF','RDFREF','POT','NEWPOT'):
        print 'Wrong DF!'
        return
    RDFsList=[]    
    name=''
    RMin=0
    RMax=0
    Npoints=0
    RDF=np.zeros((1,2))
    def _GrepText_(DFtype='inter',header='RDF ',column=2):
                '''
                Parsing set of lines <lines> for  regular expression <s> and returnd object of a proper class (it is always subclass of CDF)
                '''
                print 'Pair '+ str(pair)
                RDF_r_List=[]
                RDF_gr_List=[]
                name=basename+pair+' i'+str(irep)
                # search for resulting output (Final output, not a check one)
                for l in lines:
                  if re.search(s,l):
                    RDF_r_List.append(float(l.split()[0]))
                    RDF_gr_List.append(float(l.split()[int(column)-1]))
                RDF=np.array([[RDF_r_List[i],RDF_gr_List[i]] for i in xrange(0,len(RDF_r_List))])
                RMin=RDF[:,0].min()
                RMax=RDF[:,0].max()
                Npoints=len(RDF[:,0])
                resol=(RMax-RMin)/(Npoints-1)
                RMin=RMin-resol/2
                RMax=RMax+resol/2
                for i in xrange(0,Npoints-1):
                    if (abs(RDF[i+1,0]-RDF[i,0]-resol)>=0.0001):
                        print 'Error in RDF, points are not equidistant',RMax,RMin,Npoints,resol,RDF[0,i+1],RDF[0,i]
                        return                
                if (DFtype=='intra') and (header=='Pot1-3 '):                    
                    trdf=CDF.C_Pot13(DFtype,header+name,RMin,RMax,Npoints,g=RDF,resol=resol)
                elif (DFtype=='intra') and (header=='Pot '):
                    trdf=CDF.C_PotBond(DFtype,header+name,RMin,RMax,Npoints,g=RDF,resol=resol)
                elif (DFtype=='inter') and (header=='Pot '):
                    trdf=CDF.C_PotInter(DFtype,header+name,RMin,RMax,Npoints,g=RDF,resol=resol)
                elif (header=='RDFref ') or (header=='ADFref '):
                    trdf=CDF.C_RDFref(DFtype,header+name,RMin,RMax,Npoints,g=RDF,resol=resol)
                else:
                    trdf=CDF.C_DF(DFtype,header+name,RMin,RMax,Npoints,g=RDF,resol=resol)
                return trdf
    
    rdffile=open(rdffilename,'r')
    lines=rdffile.readlines()
    # Get total number of iterations performed:
    if (iters==()):
        TotalNumIter=0
        for l in lines:
            if re.search('Final Corrected effective potential',l):
                TotalNumIter+=1
        print 'Total number of iterations:', TotalNumIter
        iters=range(1,TotalNumIter+1)
    else:
         TotalNumIter=len(iters)
    # Get number of checks for every iteration
    TotalNumCheck=0
    for l in lines:
        if re.search('Global Potential correction convergence check',l):
            TotalNumCheck+=1
    print 'Total number of checks per iteration:', TotalNumCheck/TotalNumIter
    
    for irep in iters: #loop over all iterations
        print 'Iteration %i' %(irep)
        RDFs=[]
        for pair in PairNamesList[0]:    
            if (DF.upper()=='RDFREF')and(iters.index(irep)==0):
                s='pot: *'+pair.split('-')[0].strip()+' *'+pair.split('-')[1].strip()+' *Final *'+str(irep)+' '
                RDFs.append(_GrepText_(DFtype='inter',header='RDFref ',column=3))
            else:
            #Loop over intermolecular pairs
                s='pot: *'+pair.split('-')[0].strip()+' *'+pair.split('-')[1].strip()+' *Final *'+str(irep)+' '
                if DF.upper()=='RDF':
                    RDFs.append(_GrepText_(DFtype='inter',header='RDF '))
                if DF.upper()=='POT':
                    RDFs.append(_GrepText_(DFtype='inter',header='Pot ',column=5))
                if DF.upper()=='NEWPOT':
                    RDFs.append(_GrepText_(DFtype='inter',header='Pot ',column=4))
                if DF.upper()=='POTCORR':
                    RDFs.append(_GrepText_(DFtype='inter',header='PotCorr ',column=6))
        for pair in PairNamesList[1]:
            if (DF.upper()=='RDFREF')and(iters.index(irep)==0):
                s='pot: *'+pair.split('-')[0].strip()+' *'+pair.split('-')[1].strip()+' *'+'[0-9] *Final *'+str(irep)+' '
                RDFs.append(_GrepText_(DFtype='intra',header='RDFref ',column=3))
            else:
            #Loop over intramolecular pairs with pairwise bonds
                s='pot: *'+pair.split('-')[0].strip()+' *'+pair.split('-')[1].strip()+' *'+'[0-9] *Final *'+str(irep)+' '
                if DF.upper()=='RDF':
                    RDFs.append(_GrepText_(DFtype='intra',header='RDF '))
                if DF.upper()=='POT':
                    RDFs.append(_GrepText_(DFtype='intra',header='Pot ',column=5))
                if DF.upper()=='NEWPOT':
                    RDFs.append(_GrepText_(DFtype='intra',header='Pot ',column=4))
                if DF.upper()=='POTCORR':
                    RDFs.append(_GrepText_(DFtype='intra',header='PotCorr ',column=6))
        for pair in PairNamesList[2]:
            if (DF.upper()=='RDFREF')and(iters.index(irep)==0):
                s='pot: *'+pair.split('-')[0].strip()+' *'+pair.split('-')[1].strip()+' *'+pair.split('-')[2].strip()+' *'+'[0-9] *Final *'+str(irep)+' '
                RDFs.append(_GrepText_(DFtype='intra',header='ADFref ',column=3))
            else:
                #Loop over intramolecular triplets with 1-3 bonds
                s='pot: *'+pair.split('-')[0].strip()+' *'+pair.split('-')[1].strip()+' *'+pair.split('-')[2].strip()+' *'+'[0-9] *Final *'+str(irep)+' '
                if DF.upper()=='RDF':
                    RDFs.append(_GrepText_(DFtype='intra',header='ADF '))
                if DF.upper()=='POT':
                    RDFs.append(_GrepText_(DFtype='intra',header='Pot1-3 ',column=5))
                if DF.upper()=='NEWPOT':
                    RDFs.append(_GrepText_(DFtype='intra',header='Pot1-3 ',column=4))                    
                if DF.upper()=='POTCORR':
                    RDFs.append(_GrepText_(DFtype='intra',header='PotCorr1-3 ',column=6))
        RDFsList.append(RDFs)
    return RDFsList
def GetPotsFromFile_pot(rdffilename,basename='',mcmfile='',Ucut=10000.0):    
    '''
    Analyse pot-file <rdffilename> produced by MAGIC and return a list of potentials.
    basename - prefix given to all the potentials readed from the file
    mcmfile - optional file to read atom names, if not stated a warning will be givean, and atom numbers will be stated instead of types
    Ucut - Criteria to cut artifitial core part of the potential.
    '''
    import pdb
    PotList=[]
    rdffile=open(rdffilename,'r')
    lines=rdffile.readlines()
    line=lines.pop(0)
    NAtomTypes=int(line.split()[0])
    NPoints0=int(line.split()[1])
    NInter=int(line.split()[2])
    NIntra=int(line.split()[3])
    N1_3=int(line.split()[4])    
    RMin0=float(line.split()[5])
    RMax0=float(line.split()[6])
# Reading MCM file
    if mcmfile=='':
        print 'No mcm file given, use numbers instead of names for atoms and atomic types'
        GotTypeNames=False
        GotAtomNames=False
    else:
        mcm=open(mcmfile,'r')
        mcmlines=mcm.readlines()
        mcmline=mcmlines.pop(0).split('\n')[0]
        NAtoms=int(mcmline)
        AtomsName={}
        AtomsCharge={}
        for i in xrange(0,NAtoms):
            AtomsName[i+1]=mcmlines[i].split()[0]
            AtomsCharge[i+1]=float(mcmlines[i].split()[5])
        GotAtomNames=True
        AtomsType=[int(mcmlines[i].split()[6]) for i in xrange(0,NAtoms)] # Number of atomic type of given atom
        AtomsCharge=[float(mcmlines[i].split()[5]) for i in xrange(0,NAtoms)] # Charge of given atom
        AtomTypesCharge={} # Charge of the given atom type
        for i in xrange(0,NAtoms):
            AtomTypesCharge[AtomsType[i]]=AtomsCharge[i]
        if len(mcmlines[0].split())==8:
            GotTypeNames=True
            AtomsTypeName=[mcmlines[i].split()[7] for i in xrange(0,NAtoms)]
            AtomTypes={}        
            for i in xrange(0,NAtoms):
                AtomTypes[AtomsType[i]]=AtomsTypeName[i]
        else:
            GotTypeNames=False
# Reading Intermolecular potentials       
    if NAtomTypes>0:
        resol=(RMax0-RMin0)/NPoints0
    for iNInter in xrange(0,NInter):                
        NPoints=NPoints0
        RMin=RMin0
        RMax=RMax0        
        while float(lines[0].split('\n')[0].split()[1])>Ucut:             
            line=lines.pop(0).split('\n')[0]
            RMin+=resol            
            NPoints-=1
        RDFs=np.zeros((NPoints,2))
        for i in xrange(0,NPoints):                                                    
            line=lines.pop(0).split('\n')[0]
            r=float(line.split()[0])
            g_r=float(line.split()[1])
            ai=int(line.split()[2])
            aj=int(line.split()[3])
            ind=int((r-RMin)/resol)
            RDFs[ind,0]=r
            RDFs[ind,1]=g_r             
        #Creating CDF.PotInter - object
        if GotTypeNames:
            tpot=CDF.C_PotInter('inter',basename+'Pot '+AtomTypes[ai]+'-'+AtomTypes[aj],RMin,RMax,NPoints,g=RDFs[:,:])
        else:
            tpot=CDF.C_PotInter('inter',basename+'Pot AtomType'+str(ai)+'-AtomType'+str(aj),RMin,RMax,NPoints,g=RDFs[:,:])
        tpot.qq=AtomTypesCharge[ai]*AtomTypesCharge[aj]
        PotList.append(tpot)
# Reading Intramolecular pair potentials
    for i in xrange(0,NIntra):
        line=lines.pop(0).split('\n')[0]
        npairs=int(line.split()[0])
        ai=int(line.split()[1])
        aj=int(line.split()[2])
        NPoints=int(line.split()[3])
        RMin=float(line.split()[4])
        RMax=float(line.split()[5])
        resol=(RMax-RMin)/NPoints
        for i in xrange(0,npairs):
            line=lines.pop(0)
        RDFs=np.zeros((NPoints,2))
        for j in xrange(0,NPoints):
            line=lines.pop(0).split('\n')[0]
            r=float(line.split()[0])
            g_r=float(line.split()[1])
            ind=int((r-RMin)/resol)
            RDFs[ind,0]=r
            RDFs[ind,1]=g_r 
        #Creating CDF.PotBond - object
        if GotAtomNames:
            tpot=CDF.C_PotBond('intra',basename+'Pot '+AtomsName[ai]+'-'+AtomsName[aj],RMin,RMax,NPoints,g=RDFs[:,:])
        else:
            tpot=CDF.C_PotBond('intra',basename+'Pot Atom'+str(ai)+'-Atom'+str(aj),RMin,RMax,NPoints,g=RDFs[:,:])
        PotList.append(tpot)
# Reading Intramolecular 1-3 angle bending potentials
    for i in xrange(0,N1_3):
        line=lines.pop(0).split('\n')[0]
        npairs=int(line.split()[0])
        ai=int(line.split()[1])
        aj=int(line.split()[2])
        ak=int(line.split()[3])
        NPoints=int(line.split()[4])
        RMin=float(line.split()[5])
        RMax=float(line.split()[6])
        resol=(RMax-RMin)/NPoints
        for i in xrange(0,npairs):
            line=lines.pop(0)
        RDFs=np.zeros((NPoints,2))
        for j in xrange(0,NPoints):
            line=lines.pop(0).split('\n')[0]
            r=float(line.split()[0])
            g_r=float(line.split()[1])
            ind=int((r-RMin)/resol)
            RDFs[ind,0]=r
            RDFs[ind,1]=g_r 
        #Creating CDF.Pot13 - object
        if GotAtomNames:
            tpot=CDF.C_Pot13('intra',basename+'Pot13 '+AtomsName[ai]+'-'+AtomsName[aj]+'-'+AtomsName[ak],RMin,RMax,NPoints,g=RDFs[:,:])
        else:
            tpot=CDF.C_Pot13('intra',basename+'Pot13 Atom'+str(ai)+'-Atom'+str(aj)+'-Atom'+str(ak),RMin,RMax,NPoints,g=RDFs[:,:])
        PotList.append(tpot)
    return PotList

