import warnings
warnings.filterwarnings("ignore")
import MDAnalysis
import MDAnalysis.analysis.rdf
import numpy as np
from pymatgen.core import lattice
import builtins
import os

def rdf(traject,latt,pair,endrange,resolution):
    cell = []
    with open(latt, 'r') as f:
        for line in f:
            if 'TV' in line:
                cell.append(list(map(float, line.split()[1:4])))

    box = lattice.Lattice(cell)
    (lx, ly, lz), (a, b, c) = box.lengths_and_angles

    u = MDAnalysis.Universe(traject, format='XYZ')
    u.dimensions = np.array([lx, ly, lz, a, b, c], dtype=np.float32)

    builtins.u = u

    ag1 = None
    ag2 = None
    elems = pair.split('-')

    # Definisi file keluaran
    outfile = open('rdf_{}-{}.dat'.format(elems[0],elems[1]), 'w')

    ag1 = u.select_atoms(f'name {elems[0]}')
    ag2 = u.select_atoms(f'name {elems[1]}')

    nbins = int((endrange - 0.05)/resolution)

    rdf = MDAnalysis.analysis.rdf.InterRDF(ag1, ag2, nbins=nbins, range=(0.05, endrange),   verbose=True)
    rdf.run()

    vol = np.power(rdf.edges[1:], 3) - np.power(rdf.edges[:-1], 3)
    vol *= 4/3.0 * np.pi   
    box_vol = box.volume
    pair_den = box_vol / (len(ag1.atoms) * len(ag2.atoms))
    integral = rdf.rdf / (pair_den / vol)

    total = 0.0
    step = 0
    for r, v, d in zip(rdf.bins, rdf.rdf, integral):
        step += 1
        total += d/len(ag1.atoms)
        if (step == 1):
            continue
        print('{:8.4f} {:8.4f} {:8.4f}'.format(r, v, total), file=outfile)

    outfile.close()

    with open('rdf_{}-{}.gp'.format(elems[0],elems[1]),'w') as fout:
        print("""set terminal pdf
        set output "rdf_{}-{}.pdf"
        set key
        set yrange [0:10]
        set ytics (0,1,2,3,4,5,6,7,8,9,10)
        set xlabel "Distance [Angstroms]"
        set ylabel "RDF [arb.]"
        set style line 1 lc rgb '#00008b' lt 1 lw 2 pt 7
        set style line 2 lc rgb '#FF0000' lt 1 lw 2 pt 7
        plot 'rdf_{}-{}.dat' u 1:2 with line title "RDF" linestyle 1,\
            '' u 1:3 with line title "Coordination Number" linestyle 2""".format(elems[0],elems[1],elems[0],elems[1]),file=fout)
    os.system('gnuplot rdf_{}-{}.gp'.format(elems[0],elems[1]))