OVITO > Support Forum

Calculating RDF for molecules, not atoms

(1/2) > >>

Dear all,

I need to calculate the Radial Distribution Function for the center of the mass of molecules as a reference point. I know that I should use the modifier "Coordination Analysis", but the problem is that as I know, this modifier only considers atoms and performs the calculations. And what I want is not to the calculations atom by atom but to some specific points, for example, the center an aromatic ring or center of the mass of the molecules. I really appreciate any help from you. Thank you in advance.



Constanze Kalcher:
Dear Sahar,

what you're asking for should be achievable by using a Python script modifier, however, that will require some python coding on your part. Can you give me an example of how your input data format looks like? Then I can give some help with the script.


Dear Constanze,

Thank you for your response. Here is a lammps trajectory file which I use as an input for ovito. This file consists of all the atoms of my system and some of them together form different molecules and I need to do the rdf calculations based on molecules.

Thank you 

Constanze Kalcher:
Dear Sahar,

here's an idea how you could perform the calculation: So each atom also has a particle property "Molecule Identifier" such that you know to which molecule it belongs. In order to compute RDFs for your molecules you first need to calculate the center of mass of each molecule. This could be done by a using a python script modifier in the graphical user interface:

--- Code: ---from ovito.data import *
import numpy as np

def modify(frame, data):

    num_particles = data.particles.count
    pos = data.particles_["Position"]
    #Part 1
    #Create new particle property that contains the atomic masses
    # Edit the following list
    my_masses = [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
    mass_property = data.particles_.create_property('Mass')
    p_types = data.particles["Particle Type"]
    for atom_index in range(data.particles.count):
        with mass_property:
            mass_property[atom_index] = my_masses[ (p_types[atom_index] - 1) ]
    #Part 2
    #Calculate center of mass for every molecule and save as new particle property
    mol_IDs = data.particles_['Molecule Identifier']
    mol_com_property = data.particles_.create_property('Molecule COM', dtype = float, components = 3)
    #Loop over all Molecule Identifiers
    for current_molecule in set(mol_IDs):
        #Find all the atoms that belong to that molecule and average their positions
        with mol_com_property:
            mol_com_property[(mol_IDs == current_molecule)] = np.average( pos[((mol_IDs == current_molecule))], weights = mass_property[((mol_IDs == current_molecule))], axis = 0)
--- End code ---

This custom modifier has several parts. In part 1 you will create a particle property "Mass" and add the corresponding atomic masses since they are (currently) not automatically read in when importing a dump file. In this example I used a look up list called "my_masses" that contains the masses of all 11 particle types (they're all 1 in this test case), which you should edit.
This information is necessary to compute the center of mass vector for each molecule, which is calculated and stored in a new particle property 'Molecule COM' in Part 2 (see result in the attached screenshot).

From there you can go on and compute RDF's by either writing your own script or maybe a simpler procedure would be to export the molecule center of masses together with the molecule types (to the lammps dump or xyz format e.g.) as if they were atoms. In that way you can re-import that file and simply apply the coordination analysis modifier in the graphical user interface.

Please note, however, that the above script does not take care of periodic boundaries yet, so you'll need to add that part where you average the positions in case you have applied periodic boundary conditions in you simulation. Also, can you explain how you defined the molecule identifier? I noticed that in some cases the atoms that share the same  molecule identifier, e.g. molecule 1, are very far apart. See screenshot 2 for an example. Was that done deliberately?

Hope that helps,


Thank you very much for your response. About the molecule IDs, yes that is wrong because when I make my data files, the numbering of the molecules restart from 1 for every different component and they are not continuous for all the components of my system.


[0] Message Index

[#] Next page

Go to full version