Forum Navigation
You need to log in to create posts and topics.

Inter and Intera molecular RDF(Radial Distribution Function)

Hi,

I am trying to write an script to calculate Inter and Intera molecular RDF for polymer molecules (single type monomers). however from reading source code I assume "CutoffNeighborFinder" is only returns particle_types and not Molecule_identifier with the neighbor list particles. If it does please let me know how I can access to that data in python script. is there any ready script for my purpose? any recommendation to make the script easy?

Your response is very much appreciated.

best,

Matt

Hi Matt,

the find()function of the utility class ovito.data.CutoffNeighborFinder returns an iterator over all neighbors of the given particle. Each of these neighbors is an object with attributes, such as e.g. its particle index index and the distance of the current particle to the central distance. Once you know the index of a particle you can look up any of its properties. That means you could in principle use the particle index of each neighbor to look up the Molecule ID of that neighbor particle.

You could use the code example from the documentation, https://www.ovito.org/docs/current/python/modules/ovito_data.php#ovito.data.CutoffNeighborFinder,
and substitute the particle type property with the Molecule ID property. Assuming the name of the particle property is "Molecule Identifier", then you can look up the Molecule ID of each neighbor particle like this:

mtypes = data.particles["Molecule Identifier"]
# Loop over all particles:
for index in range(data.particles.count):
    # Iterate over the neighbors of the current particle:
    for neigh in finder.find(index):
        # The index can be used to access properties of the current neighbor, e.g.
        molecule_id_of_neighbor_particle = mtypes[neigh.index]
        print(neigh.index, molecule_id_of_neighbor_particle)

I'm wondering though if it could be easier to use the CoordinationAnalysisModifier? If I'm not misunderstanding the problem setup, you could use this out-of-the-box solution to generate the intra-molecular RDF's. To calculate inter-molecular RDF's, a trick could be to calculate the center of mass (COM) of each molecule and then export the molecule id and molecule type along with the COM-coordinates as "generalized" particle objects so to say. That way you can again apply the CoordinationAnalysisModifier to that data without having to implement too much extra functionality yourself.

Let me know what you think. I'm also happy to help with developing the python script.

-Constanze

Dear Constanze,

Thank you for your helpful response.

I thought about using CoordinationAnalysisModifier, and what I have now so far is attached. but I guess I cannot use it since it calls NeighborFinder itself again and resets the neighbor lists.

I guess I have to forget about changing the main neighbor lists data and just add an script to calculate g(r).

Also about your suggestion I guess replacing a molecule COM with it-self and calculation of intra-G(r) only would work for small molecules.

Thank you,

M.

 

Uploaded files:
  • You need to login to have access to uploads.

Hi,

So I think I have something (attached file) that might be what I need. but I am getting some errors.

First, why "import_file" does not read lammps dump file?

second, is there any "ovito.pipeline"?  as in https://www.ovito.org/docs/current/python/modules/ovito_pipeline.php#ovito.pipeline.Pipeline.source  I have the error {ImportError: No module named 'ovito.pipeline'}

how I can use the data that is already loaded in Ovito?

I would appreciate your help to solve the errors. although, I am not sure about the normalization part but it is something that can be fixed later if script works.

Thank you,

Matt

 

Uploaded files:
  • You need to login to have access to uploads.

Hi Matt,

could you maybe share an example of your input data with me (either here or through the official mail support)? That would make it easier for me to help you.
The import_file statement in your script looks correct, so I can only guess it's either an issue with the file path or the dump file?
Another situation in which this error message could arise that you were maybe trying to use this script in a python script modifier in the GUI? Please note that you have to execute your script from the terminal by calling ovitos and that it can't be used in a Python script modifier in the GUI. If that was the case please have a look at this introductory section in the scripting manual: https://www.ovito.org/docs/current/python/introduction/running.php
You might also find this article helpful that explains OVITO's data model: https://www.ovito.org/docs/current/python/introduction/data_model.php#data-model-intro,
therein, you'll find a section on how to look up particle properties. That was what your question "how I can use the data that is already loaded in Ovito" was aiming at, was it?

One thing I do not understand from your script, though, is this line:
g_r, bin_edges=np.histogram(finder.find(distance),bin_num).
Firstly, you need to pass a particle index to the finder.find() and secondly, there is no variable distance defined in your script.

Here is my suggestion for generating a histogram of all particle distances, i.e. of particle pairs that do not belong to the same molecule.

num_bins = 100
bin_range = (0, cutoff)
g_r = numpy.zeros(num_bins)
mtypes = data.particles["Molecule Identifier"]
# Loop over all particles:
for index in range(data.particles.count):
    molecule_id_of_center_particle = mtypes[index]
    # Iterate over the neighbors of the current particle:
    for neigh in finder.find(index):
        # The index can be used to access properties of the current neighbor, e.g.
        molecule_id_of_neighbor_particle = mtypes[neigh.index]
        if molecule_id_of_neighbor_particle != molecule_id_of_center_particle:
            g_r += numpy.histogram(neigh.distance, range = bin_range, bins = num_bins )[0]
#Normalize g_r
#...

#Export
r = numpy.linspace(0, cutoff, num=num_bins, endpoint=False)
numpy.savetxt("my_rdf.txt", numpy.column_stack((r , g_r)))

Let me know if you would like to work in the GUI instead, since it is possible to plot and display your computation results in the Data tables tab of the Data inspector. https://www.ovito.org/docs/current/python/modules/ovito_data.php#ovito.data.DataTable

-Constanze

Constanze,

Thank you for your help. I've completed the code. I will upload the code here so people can use it, As soon as I'm sure the normalization is correct.

Best Wishes and Be healthy

MB

Great! I'm happy to hear that you could make good progress with this.