Author Topic: how to count the number of B type neighbors to certain A type atoms  (Read 71 times)

kozak

  • Newbie
  • *
  • Posts: 2
Hi there,
I am a newbie to Ovito, please help me if you can! ;)
I met a probloem, I am using a polymer system with 30000+ atoms with two types (A: carbon atom, B: hydrogen atom). Now I am only intersted in one carbon atom with certain ParticleIdentifier (say 21460), I want to count the number of neighbors (only hydrogen atoms) to this carbon atom . In python script, carbon atoms are type one and hydrogen atoms are type two, the cutoff distance to build neighbors between hydrogen and carbon atoms is 1.30. In this case, if the coordination number of hydrogen atom is equal to one, which means it bonded to the carbon atom, and I want to count the number of this hydrogen atoms. The results should be 1 or 2 or 3 maybe. With some previous references, I wrote a script but it didnt work. I think the reason is that the modifiers I used after selecting the carbon atom always calculate all atoms in the system. But in my case I want it to perform in one carbon atom and all hydogen atoms.

import ovito
from ovito.io import import_file
from ovito.modifiers import *
import numpy as np

# Load in the data file
node = import_file("1.dump", multiple_frames = True)

# Output the number of atoms and frames
print("Number of atoms = %i" % node.source.number_of_particles)
print("Number of frames = %i" % node.source.num_frames)

node.modifiers.append(SelectExpressionModifier(expression = 'ParticleIdentifier == 18267'))
node.compute()

# Calculate the coordination number of the hydrocarbon system, build neighbour C-H bonds
modifier = ComputePropertyModifier()
modifier.cutoff_radius = 1.30
modifier.output_property = "Coordination"
modifier.expressions = ["0"]
modifier.neighbor_mode = True
modifier.neighbor_expressions = ["ParticleType==2?1:0"]
node.modifiers.append(modifier)

# Add a user-defined modifier that counts the number of particles of type 2
# and which have a coordination equal to 1.
def count_coordinated_particles(frame, input, output):
   n_hydrogen = (input.particle_properties['Particle Type'].array == 2) & (input.particle_properties['Coordination'].array == 1)
   output.attributes['coordinated_particles'] = np.count_nonzero(n_hydrogen)

# Insert Python modifier into the data pipeline.
node.modifiers.append(PythonScriptModifier(function = count_coordinated_particles))

# Calculate the neighbours for the last frame only
node.compute(node.source.num_frames)
   
print(node.output.attributes['coordinated_particles'])

Again, please help me with it and thank you so so much!!! :)

Constanze Kalcher

  • Administrator
  • Newbie
  • *****
  • Posts: 35
Re: how to count the number of B type neighbors to certain A type atoms
« Reply #1 on: August 01, 2018, 05:55:20 PM »
Hi,

actually, here you calculate the number of neighbors of type 2 within a cutoff of 1.3 for every atom.
Quote
Code: [Select]
modifier = ComputePropertyModifier()
modifier.cutoff_radius = 1.30
modifier.output_property = "Coordination"
modifier.expressions = ["0"]
modifier.neighbor_mode = True
modifier.neighbor_expressions = ["ParticleType==2?1:0"]
node.modifiers.append(modifier)

If you only need that information for the particle you previously selected you also need to use
Code: [Select]
only_selected = True in your Compute Property modifier.
In that case, the particles will have a particle property called "Coordination" which only for the selected particle should be equal to the number of type2 neighbors (within that cutoff) and zero for all other particles.

Then in the following part of your script you count the number of all hydrogen atoms that have exactly one hydrogen neighbor. What were you trying to achieve here?
Quote
Code: [Select]
def count_coordinated_particles(frame, input, output):
   n_hydrogen = (input.particle_properties['Particle Type'].array == 2) & (input.particle_properties['Coordination'].array == 1)
   output.attributes['coordinated_particles'] = np.count_nonzero(n_hydrogen)

-Constanze
« Last Edit: August 01, 2018, 06:19:07 PM by Constanze Kalcher »

kozak

  • Newbie
  • *
  • Posts: 2
Re: how to count the number of B type neighbors to certain A type atoms
« Reply #2 on: August 02, 2018, 01:02:29 AM »
Thank you for the quick reply!
I want to know how many hydrogen atoms bonded to the selected carbon atom. (This is used to determine the free radiacal) 
Can you please write me a syntax that can do the statistics of that? That would be so great! Thank you for your help!

Constanze Kalcher

  • Administrator
  • Newbie
  • *****
  • Posts: 35
Re: how to count the number of B type neighbors to certain A type atoms
« Reply #3 on: August 03, 2018, 10:47:37 AM »
Hi,
in principle you were only missing a loop over all your frames, which I added in the end. Also I removed a couple of node.compute()'s from your code. You don't need to apply it after every modifier addition. You can just do it after you set up your modification pipeline. :)

Try the following code (it's still OVITO 2 syntax) which calculates the average hydrogen coordination number of your selected particle and let me know if you have questions:

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

#Step 1: Setting up the pipeline
# Load in the data file
node = import_file("1.dump", multiple_frames = True)

# Output the number of atoms and frames
print("Number of atoms = %i" % node.source.number_of_particles)
print("Number of frames = %i" % node.source.num_frames)

node.modifiers.append(SelectExpressionModifier(expression = 'ParticleIdentifier == 18267'))

# Calculate the coordination number of the hydrocarbon system, build neighbour C-H bonds
modifier = ComputePropertyModifier()
modifier.cutoff_radius = 1.30
modifier.output_property = "Coordination"
modifier.expressions = ["0"]
modifier.neighbor_mode = True
modifier.neighbor_expressions = ["ParticleType==2?1:0"]
modifier.only_selected=True
node.modifiers.append(modifier)

n_hydrogen = []
# Step 2: Evaluating the same pipeline for different input frames:                                                                                                                                                                 
for frame in range(node.source.num_frames):
        data = node.compute(frame)
        coordination = data.particle_properties["Coordination"].array
        selection = data.particle_properties["Selection"].array
        n_hydrogen.append( coordination[ (selection == 1) ]  )

print("average hydrogen coordination number of selected particle: %.2f" % np.mean(n_hydrogen))

-Constanze