The link to the trajectory is here.I am using the module and similar code of the ASAP calculator from ASE. The problem is as follows, the running Coordination Number doesn't match the number of atoms in the cell. Upon checking, I found that the ASE is reading the trajectory atoms boundary condition and cell correctly. But the results are this with 16 Ge-atoms in the cell. enter image description here enter image description here enter image description here

So, I am curious if someone can get the simple code below working,

from scipy import integrate
from asap3 import *
from asap3.analysis.rdf import RadialDistributionFunction
from ase.lattice.compounds import *
import matplotlib.pyplot as plt
import numpy as np
from ase.data import atomic_numbers, atomic_names, atomic_masses, covalent_radii

def rdf(file_traj, ge_ratio):
    rMax = 12.41
    nBins = 200
    traj = Trajectory(file_traj)
    RDFobj = None
    for atoms in traj:
        if RDFobj is None:
            RDFobj = RadialDistributionFunction(atoms, rMax, nBins)
            RDFobj.atoms = atoms  # Fool RDFobj to use the new atoms
        RDFobj.update()           # Collect data

    x = (np.arange(nBins) + 0.5) * rMax / nBins      
    rdf = RDFobj.get_rdf()       
    plt.title('Total RDF-Ge Count%s'%ge_ratio, fontsize=15)
    plt.plot(x, rdf)

   # Get the partial RDFs and plot them
    Si = atomic_numbers['Si']
    Ge = atomic_numbers['Ge']

    rdfGeGe = RDFobj.get_rdf(elements=(Ge, Ge))
    plt.title('Partial RDF-Ge Count%s'%ge_ratio, fontsize=15)
    plt.plot(x, rdfGeGe, 'r', label = "Ge-Ge")
    y = rdfGeGe * 4 * np.pi * x * x * (ge_ratio/rMax**3)   
    result_lst = []   
    for i in range(len(x)):
        result = integrate.simps(y[:i + 1],x[: i+ 1])
    result_lst = np.array(result_lst)

    plt.plot(x, result_lst, 'r', label = "Ge-Ge")
    plt.title('CN-Ge Count%s'%ge_ratio, fontsize=15)
  • $\begingroup$ This is from their documentation about the method:--- IMPORTANT: The partial RDFs are normalized such that they sum up to the global RDF. This means that integrating the first peak of a partial RDF obtained with elements=(a,b) does not give you the number of B atoms in the first shell around the A atom. Instead it gives this coordination number multiplied with the concentration of A atoms $\endgroup$
    – Saha_1994
    Dec 14, 2021 at 7:53
  • $\begingroup$ Can we have the trajectory file that you are reading in? This will allow users to try to correct your code using your actual data to test it. $\endgroup$
    – Tyberius
    Dec 14, 2021 at 16:00
  • $\begingroup$ The Google drive link for the trajectory file, drive.google.com/file/d/1TwkQQADbhKsS59wxgVte3Hcz_3KugSpK/… Also, to use the func rdf, use rdf(".YOURPATH/local_minima.traj", 16) $\endgroup$
    – Saha_1994
    Dec 14, 2021 at 17:11
  • $\begingroup$ @Saha_1994 I'm unable to view that file in Google Drive, can you please upload it to a folder here called 7256 (since that's this question's ID number in its URL)? $\endgroup$ Dec 15, 2021 at 7:52
  • 1
    $\begingroup$ @NikeDattani While it would probably be better to have on GitHub, I think the no preview is just because these .traj files are binary and Google Drive doesn't have a default way to display a generic binary file. $\endgroup$
    – Tyberius
    Dec 20, 2021 at 15:58