Author Topic: Distribution of coordination number  (Read 1873 times)

Mahsa

  • Newbie
  • *
  • Posts: 9
Distribution of coordination number
« on: January 31, 2019, 06:25:15 PM »
Dear supports,

I am trying to get the distribution of coordination number for special atoms around Li ion in my system. I could get the average value of RDF and CN from Gromacs tools before. I converted the trajectory file from MD simulation to a pdb file and then try to analyse that in OVITO.
I followed this post in OVITO forum:
http://forum.ovito.org/index.php?topic=340.0
I did the same steps as in the post, however, the histogram does not show anything and in the data inspector the selection for all particle is 0 (Please see attached to this post the screenshots). What did I miss in the analysis?

Another question is, I have different type of O in the molecules. I'd like to calculate the distribution of coordination numbers of different types of O around Li ions, separately over the trajectory. How can I do this in OVITO?

Best regards,
Mahsa


Constanze Kalcher

  • Administrator
  • Sr. Member
  • *****
  • Posts: 299
Re: Distribution of coordination number
« Reply #1 on: February 01, 2019, 10:13:26 AM »
Dear Mahsa,

could you upload an example of your structure here and also paste the python script, then I can help you more easily.

-Constanze

Mahsa

  • Newbie
  • *
  • Posts: 9
Re: Distribution of coordination number
« Reply #2 on: February 01, 2019, 10:52:51 AM »
Dear Constanze,

I can't upload the file because the size is above the limit. Is there another way that I can share the file with you?
This is the script I used from the other post in the forum:

from ovito.data import *
import numpy as np

def modify(frame, input, output):
    # Prefetch the property array containing the particle type information:
    ptypes = input.particles['Particle Type']
   
    #Create a new particle property that stores the number of neighbors with particle type 2 and initialize with zeros
    values = np.zeros(input.particles.count, dtype=int)
    new_property = output.particles.create_property('Type2 Coordination', data=values)
   
    # Initialize neighbor finder object:
    cutoff = 3.
    finder = CutoffNeighborFinder(cutoff, input)
   
    # Loop over all particles:
    for index in range(input.particles.count):
        #Only iterate over particles of type 1
        if(ptypes[index] == 1):
           # Get a list of particle types of all neighbors of the current particle:
            neighbor_types = [ptypes[neigh.index] for neigh in finder.find(index)]
            #Count the number of neighbor types of all types and store the count for type 2 in the newly created particle property
            with new_property:
                new_property[index] = np.bincount(neighbor_types, minlength=4)[2]


-Mahsa

Constanze Kalcher

  • Administrator
  • Sr. Member
  • *****
  • Posts: 299
Re: Distribution of coordination number
« Reply #3 on: February 04, 2019, 12:03:29 PM »
Dear Mahsa,

you don't need to upload the whole trajectory, a single shapshot would be enough for me as an example. But yes you can share a download link here or sent a private message to us via the mail support.
The script doesn't seem to be the problem, are you sure the cutoff used for the CutoffNeighborFinder function is appropriate? I'm happy to look into it further once I have access to the data.

-Constanze
« Last Edit: February 04, 2019, 12:05:51 PM by Constanze Kalcher »

Mahsa

  • Newbie
  • *
  • Posts: 9
Re: Distribution of coordination number
« Reply #4 on: February 04, 2019, 12:25:54 PM »
Dear Constanze,

Please find attached to the post, 1ps of the trajectory and the converted file to pdb as an example.

The more reasonable cutoff, is 2 Å for this system but I wanted to see if the script works fine for the system.

-Mahsa

Constanze Kalcher

  • Administrator
  • Sr. Member
  • *****
  • Posts: 299
Re: Distribution of coordination number
« Reply #5 on: February 04, 2019, 02:23:48 PM »
Dear Mahsa,

alright, I had a look at your data. You were asking how to calculate the number of O-atoms (which have a particle type id = 1 as shown in the screenshot) around the Li-atoms (in your case named "Li+" and having a particle type id = 4).

In the python script you posted here you're obviously counting C atoms around O-atoms, because you iterate over particles having a particle type 1 and count how many neighbors of type 2 they have. Bottom line is, you should never just copy and paste a custom python modifier
function without understanding every bit of it.
I know it's quite a steep learning curve in the beginning, so here's the corrected version (for the latest OVITO developer version) as a head start for you:

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

def modify(frame, data):
    # Prefetch the property array containing the particle type information:
    ptypes = data.particles['Particle Type']
    #Create a new particle property that stores the number of neighbors with particle type 2 and initialize with zeros
    values = np.zeros(data.particles.count, dtype=int)
    new_property = data.particles_.create_property('Li-O Coordination', data=values)
    # Initialize neighbor finder object:
    cutoff = 2.
    finder = CutoffNeighborFinder(cutoff, data)
   
    # Loop over all particles:
    for index in range(data.particles.count):
        #Only iterate over particles of type Li
        Li_type_id = ptypes.type_by_name('Li+').id
        if(ptypes[index] == Li_type_id):
           # Get a list of particle types of all neighbors of the current particle:
            neighbor_types = np.array([ptypes[neigh.index] for neigh in finder.find(index)])
            #Count the number of  O neighbor types and store them in the newly created particle property
            O_type_id = ptypes.type_by_name('O').id
            with new_property:
                new_property[index] = np.count_nonzero(neighbor_types == O_type_id)

In case you are unsure about what particle type id corresponds to what particle type name  the function type_by_name().id is very useful.
https://ovito.org/manual_testing/python/modules/ovito_data.html#ovito.data.Property

Furthermore, you asked if you could distinguish between two different groups of O-atoms. Yes, in your case that's possible because they have a distinguishable feature which is the Molecule Type.
So you have O atoms which have a Particle Type 1 and a Molecule Type 1 ("LIG") and another group of O-atoms which have a Particle Type 1 and a Molecule Type 3 ("tf2") and you
can select the atoms that fulfill both conditions, e.g. just

add a line at the top where you get the molecule types
Code: [Select]
mtypes = data.particles['Molecule Type']
and then e.g. replace the if condition by:

Code: [Select]
if(ptypes[index] == Li_type_id):
           # Get a list of particle types of all neighbors of the current particle:
            neighbor_types = [ ( ptypes[neigh.index], mtypes[neigh.index]) for neigh in finder.find(index)]
            #Count the number of  O neighbor types and store them in the newly created particle property
            print(neighbor_types)
            O_type_id = ptypes.type_by_name('O').id
            tf2_type_id = mtypes.type_by_name('tf2').id
            print(tf2_type_id)
            with new_property:
                new_property[index] = neighbor_types.count( (O_type_id, tf2_type_id))

where you use both the Particle Type and the Molecule Type as criterion in the count function.

Let me know if you have questions,

-Constanze
« Last Edit: February 04, 2019, 02:25:49 PM by Constanze Kalcher »

Mahsa

  • Newbie
  • *
  • Posts: 9
Re: Distribution of coordination number
« Reply #6 on: February 07, 2019, 03:15:55 PM »
Dear Constanze,

Thank you very much for your help and explanation!

-Mahsa