OVITO Forum

OVITO => Frequently Asked Questions (FAQ) => Topic started by: Alexander Stukowski on May 03, 2017, 04:37:26 PM

Title: How can the partial radial distribution functions be calculated with Ovito?
Post by: Alexander Stukowski on May 03, 2017, 04:37:26 PM
OVITO's Coordination Analysis (http://ovito.org/manual/particles.modifiers.coordination_analysis.html) modifier can calculate the radial distribution function (RDF) (https://en.wikipedia.org/wiki/Radial_distribution_function) 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 (http://ovito.org/manual/particles.modifiers.coordination_analysis.html) 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):


Code: [Select]
import numpy as np

num_bins = 100

def 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.
Title: Re: How can the partial radial distribution functions be calculated with Ovito?
Post by: Alexander Stukowski on May 03, 2017, 04:41:20 PM
The steps described above can also be automated. The following example Python script loads a simulation file, computes the partial RDF for pairs of particles of type 1 and 1, and writes it to a text file. The script can be run with the ovitos (http://ovito.org/manual/python/introduction/running.html#ovito-s-python-interpreter) interpreter.

Code: [Select]
from ovito.modifiers import *
from ovito.io import *
import numpy as np

num_bins = 100
cutoff = 5.0

node = 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.volume
factor = 4./3. * np.pi * rho * output.number_of_particles

radii = 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)