Regarding rdf from HOOMD trajectory (gsd)

[Arpita Srivastava]

Dear Freud Users,

I am new to Freud and trying to compute rdf between two particles, A and A1, which are present in my system. However, I am getting a blank plot, i.e., no values at x and y axes. My script is attached below, please help me in fixing this problem.

import freud
import pickle
import gsd.hoomd
import matplotlib.pyplot as plt
import numpy as np

def snap_molecule_indices(snap):
    system = freud.AABBQuery.from_system(snap)
    num_query_points = num_points = snap.particles.N
    query_point_indices = snap.bonds.group[:, 0]
    point_indices = snap.bonds.group[:, 1]
    distances = system.box.compute_distances(
        system.points[query_point_indices], system.points[point_indices]
    )
    nlist = freud.NeighborList.from_arrays(
        num_query_points, num_points, query_point_indices, point_indices, distances
    )
    cluster = freud.cluster.Cluster()
    cluster.compute(system=system, neighbors=nlist)
    return cluster.cluster_idx


def intermolecular_rdf(
    gsdfile = "s2.gsd",
    A_name = "A",
    B_name = "A",
    start=0,
    stop=None,
    r_max=10,
    r_min=0,
    bins=100,
    exclude_bonded=True,
):
       with gsd.hoomd.open(gsdfile) as trajectory:
        snap = trajectory[0]

        if r_max is None:
            # Use a value just less than half the maximum box length.
            r_max = np.nextafter(
                np.max(snap.configuration.box[:3]) * 0.5, 0, dtype=np.float32
            )

        rdf = freud.density.RDF(bins=bins, r_max=r_max, r_min=r_min)
        type_A = snap.particles.typeid == snap.particles.types.index(A_name)
        type_B = snap.particles.typeid == snap.particles.types.index(B_name)

        if exclude_bonded:
            molecules = snap_molecule_indices(snap)
            molecules_A = molecules[type_A]
            molecules_B = molecules[type_B]

        for snap in trajectory[start:stop]:

            A_pos = snap.particles.position[type_A]
            if A_name == B_name:
                B_pos = A_pos
                exclude_ii = True
            else:
                B_pos = snap.particles.position[type_B]
                exclude_ii = False

            box = snap.configuration.box
            system = (box, A_pos)
            aq = freud.locality.AABBQuery.from_system(system)
            nlist = aq.query(
                B_pos, {"r_max": r_max, "exclude_ii": exclude_ii}
            ).toNeighborList()

            if exclude_bonded:
                pre_filter = len(nlist)
                indices_A = molecules_A[nlist.point_indices]
                indices_B = molecules_B[nlist.query_point_indices]
                nlist.filter(indices_A != indices_B)
                post_filter = len(nlist)

            rdf.compute(aq, neighbors=nlist, reset=False)
#            normalization = post_filter / pre_filter if exclude_bonded else 1
#            rdf, normalization = intermolecular_rdf(gsdfile, A_name, B_name, r_max=6, exclude_bonded=False)
#            print(f"Normalization for bonded pair exclusion: {normalization:.5f}")
            return rdf
        plt.plot(rdf.bin_centers, rdf.rdf)
plt.xlabel("r (A.U.)")
plt.ylabel("g(r)")
plt.show()

Any suggestions will be of great importance to me. Thanks in advance!

-Arpita Srivastava

[Tommy Waltmann]

Hi Arpita,

Your `intermolecular_rdf` function returns before it gets to the `plt.plot(...)` statement. Also, there is no call to the `intermolecular_rdf` function, so it never executes in the first place. When you run `python script.py` with your script, the only 3 lines that get executed are

```

plt.xlabel(...)

plt.ylabel(...)

plt.show()

```

Best,

Tommy

[Arpita Srivastava]

Hi Tommy,

Thank you so much for your kind help.

Regards,

Arpita

[Arpita Srivastava]

Hi,

I have again tried to modify my code (given above), but somehow my hoomd-blue trajectory file is (s2.gsd) is not being read by the program. I am seeking guidance if I am missing something in my code or providing inputs in an inappropriate manner? I really appreciate your suggestions for helping me to modify my script and making it work well.

Thanks,

Arpita

[Tommy Waltmann]

Arpita,

I do not know how you have modified the code you posted in your original message, so I can't really help you unless you post your new script here. It would also be helpful if you could attach your gsd file, so I can reproduce the issue you are having.

Best,

Tommy

[Arpita Srivastava]

Hi Tommy,

I hereby attach my gsd file and script too where I am trying to compute rdf between atoms A and A1. Please suggest me if I am doing something incorrect within my script.

Thank you.

[Arpita Srivastava]

The gsd file was unsent in the previous email. Kindly find it here.

Thanks,

-Arpita

[Tommy Waltmann]

Arpita,

Your python script has an `intermolecular_rdf` function defined, but there is no code calling it, so it will not execute when you run your script. If you want your code to run the `intermolecular_rdf` function, you should add some lines at the end of your script that read

```

if __name__ == "__main__":

    intermolecular_rdf()

```

-Tommy