### Author Topic: How can the partial radial distribution functions be calculated with Ovito?  (Read 4151 times)

#### Alexander Stukowski

• Sr. Member
• Posts: 435
##### How can the partial radial distribution functions be calculated with Ovito?
« on: May 03, 2017, 04:37:26 PM »
OVITO's Coordination Analysis modifier can calculate the radial distribution function (RDF) of a system of particles. The RDF g(r) measures the probability of finding an atom at distance r given that there is an atom at position 0; it is essentially a normalized histogram of interatomic distances -and is calculated as such.

However, the Coordination Analysis modifier only calculates the global RDF ignoring the types of particles. For systems with more than one species,  the so-called partial radial distribution functions gαβ(r) may be of interest. To obtain them, a different route is currently needed (until the Coordination Analysis modifier will be extended in a future program version):

• Insert the Create Bonds modifier. Let it create bonds only between those types of particles for which the partial RDF should be calculated. The cutoff parameter defines the range over which the RDF will be computed.
• Insert the Compute bond lengths modifier.
• Insert a Python script modifier and copy/paste the following script function, which computes the RDF from the bond lengths:

Code: [Select]
`import numpy as npnum_bins = 100def modify(frame, input, output): hist, bin_edges = np.histogram(input.bond_properties.length.array, bins=num_bins) rho = input.number_of_particles / input.cell.volume factor = 4./3. * np.pi * rho * input.number_of_particles radii = bin_edges[:-1] radii_right = bin_edges[1:] rdf = hist / (factor * (radii_right**3 - radii**3)) print(np.column_stack((radii,rdf)))`
Note that the partial RDF calculated by this script is normalized such that all partial RDFs, when computed for all possible pairs of atom types, sum up to the global RDF.

#### Alexander Stukowski

`from ovito.modifiers import *from ovito.io import *import numpy as npnum_bins = 100cutoff = 5.0node = import_file("simulation.0.dump")node.add_to_scene()create_bonds_modifier = CreateBondsModifier(cutoff=cutoff, mode=CreateBondsModifier.Mode.Pairwise)create_bonds_modifier.set_pairwise_cutoff('Type 1', 'Type 1', cutoff)node.modifiers.append(create_bonds_modifier)node.modifiers.append(ComputeBondLengthsModifier())output = node.compute()hist, bin_edges = np.histogram(output.bond_properties.length.array, bins=num_bins)rho = output.number_of_particles / output.cell.volumefactor = 4./3. * np.pi * rho * output.number_of_particlesradii = bin_edges[:-1]radii_right = bin_edges[1:]rdf = hist / (factor * (radii_right**3 - radii**3))result = np.column_stack((radii,rdf))print(result)np.savetxt("partial_rdf.dat", result)`