Author Topic: WS Analysis on system that is moving  (Read 1665 times)

jhart

  • Newbie
  • *
  • Posts: 22
WS Analysis on system that is moving
« on: September 01, 2017, 05:54:33 AM »
Is there a way to perform ws analysis on a system who's center of mass is changing slightly? I want to just compare atoms that have moved lattice site, but WS analysis does not work once the atoms have shifted.

Alexander Stukowski

  • Administrator
  • Hero Member
  • *****
  • Posts: 638
Re: WS Analysis on system that is moving
« Reply #1 on: September 01, 2017, 12:05:30 PM »
OVITO doesn't have a built-in feature for eliminating such a drift of the center of mass. But I assume it is not too difficult to write a Python script modifier which does that.

The user-defined modifier function would first have to compute the current center of mass and determine its displacement from the reference position. In a second step, the modifier function would subtract that displacement vector from all atomic positions.

This Python script modifier has to precede the WS analysis modifier in the data pipeline so that the WS modifier sees the corrected atom coordinates.

jhart

  • Newbie
  • *
  • Posts: 22
Re: WS Analysis on system that is moving
« Reply #2 on: September 01, 2017, 06:00:27 PM »
Any tips on how to calculate the COM with python? Is there a way to just create new atomic positions by just adding x,y,z to the old positions?
« Last Edit: September 01, 2017, 06:19:00 PM by jhart »

Alexander Stukowski

  • Administrator
  • Hero Member
  • *****
  • Posts: 638
Re: WS Analysis on system that is moving
« Reply #3 on: September 01, 2017, 06:24:00 PM »
Yeah, that is the tricky part. Some atomic coordinates might wrap around when they cross the periodic boundaries of the box. This makes the calculation of the COM difficult. You are using PBCs, right?

A strategy to deal with this could be to use only a subset of atoms and exclude atoms from the COM calculation that might cross a boundary.

Alexander Stukowski

  • Administrator
  • Hero Member
  • *****
  • Posts: 638
Re: WS Analysis on system that is moving
« Reply #4 on: September 01, 2017, 06:30:30 PM »
Here is a sketch for a Python modifier function:

Code: [Select]
from ovito.data import ParticleProperty

def modify(frame, input, output):
    old_pos_property = input.particle_properties['Position']
    old_pos = old_pos_property.array

    # Somehow compute 'shift_vector' from the current positions...

    # Compute new positions by subtracting the same shift vector from all current positions:
    new_pos = old_pos - shift_vector

    # Output the coordinates as 'Position' particle property, effectively replacing old property:
    output.create_particle_property(ParticleProperty.Type.Position, data=new_pos)       
« Last Edit: September 01, 2017, 06:34:34 PM by Alexander Stukowski »

jhart

  • Newbie
  • *
  • Posts: 22
Re: WS Analysis on system that is moving
« Reply #5 on: September 18, 2017, 09:05:47 PM »
I should first begin my explaining fully what is happening. I am running a radiation damage cascade for several hundred ps. For the first 150 ps, the WS behaves as I would expect. The number of defects spikes at around 0.2 ps and then anneals. But once the simulation reaches ~150 ps, I get defects all over the simulation cell that are not even close to the cascade itself. I am assuming it is because the COM is shifting by I can't be sure. Have you ever seen this happen before? I cannot seem to figure out what is causing this. When I zoom in and look at the atoms that are supposedly vacancies, I see that they are not, just the WS thinks they are, which made me think it is because the COM has shifted. Also, note that my system is not "exploding" with energy or anything like that, the average temperature is 900 K and has not changed.

In order to remedy all of this, I have been trying to compute the shift vector as you said above in your code example.

However, in order to start this process, instead of calculating the shift vector which would be more involved, I just decided to "Create" the shift vector. This is not physical but I just want to see if it affects the results. However, when I try something like:

shift_vector=[0,15,0]

my results are identical. The shift vector is indeed changing the y-pos of all the atoms by 15 Angstroms correct?

My full code is here:

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

node = import_file("cascade.dump",multiple_frames = True)

# Perform Wigner-Seitz analysis:
ws = WignerSeitzAnalysisModifier(
    per_type_occupancies = True)
ws.reference.load("../position.dump")
node.modifiers.append(ws)

# Define a modifier function that selects sites of type A=1 which
# are occupied by exactly one atom of type B=2.
def modify(frame, input, output):
    old_pos_property=input.particle_properties['Position']
    old_pos=old_pos_property.array
    shift_vector=[0,15,0]
    new_pos=old_pos-shift_vector
    output.create_particle_property(ParticleProperty.Type.Position, data=new_pos)


    # Retrieve the two-dimensional Numpy array with the site occupancy numbers.
    occupancies = input.particle_properties['Occupancy'].array
   
    # Get the site types as additional input:
    site_type = input.particle_properties.particle_type.array

    # Calculate total occupancy of every site:
    total_occupancy = np.sum(occupancies, axis=1)

    # Set up a particle selection by creating the Selection property:

    selection1 = (site_type == 1) & (occupancies[:,0] == 0) & (occupancies[:,1] == 0)
   
    output.attributes['O_Vac'] = np.count_nonzero(selection1)

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

# Let OVITO do the computation and export the number of identified
# antisites as a function of simulation time to a text file:
export_file(node, "defects_shift.txt", "txt",
    columns = ['Timestep', 'O_Vac'],
    multiple_frames = True)
[/code]




Also, could let me know how the WS analysis is defined and what is defined as a vacancy? i.e. how far must the closest atom be from a lattice site to be considered a vacancy? If I could play with that parameter it may solve my problems here. The reason I say this is that the atoms that become Frenkel pairs randomly are necessarily at the boundaries.

Any help would be greatly appreciated!



« Last Edit: September 19, 2017, 01:03:45 AM by jhart »

Alexander Stukowski

  • Administrator
  • Hero Member
  • *****
  • Posts: 638
Re: WS Analysis on system that is moving
« Reply #6 on: September 26, 2017, 05:53:38 PM »
Yes, the high number of defects you describe sounds like an artefact of the WS analysis.

It is easy to demonstrate, as an example, that the WS analysis reports a large number of Frenkel pair defects for a perfect crystal if you rigidly shift that crystal lattice by some amount. Once the atoms lie halfway between two reference lattice sites, then they will be randomly assigned to either of the sites. Then you get a lot of sites with zero occupancy and interstitial sites.

According to the WS analysis, a vacancy is a site with no physical atom being closer to that site than to any of the other sites. One can describe this method also in terms of the Voronoi tessellation. Each site has a Voronoi cell around it, with the Voronoi faces being located halfway between the site and its neighbors. The WS analysis counts how many physical atoms are located in each Voronoi cell. Those sites whose Voronoi cell contains no atoms at all, are called vacancies.
Note that there is no threshold parameter that would control how far an atom can move away from its site before it is counted as a defect. The WS analysis is a parameter-free method. An atom may be counted as a defect as soon as it gets closer to a neighbouring site than to its original site.

So, for the WS analysis to work correctly, it is important that perfect crystal atoms are still remain close to the center of the Voronoi cell (i.e. close to the original site locations). If all the atoms somehow get shifted toward the borders of their Voronoi cells, then they might spuriously be assigned to the cell of a neighbouring site, especially if they are also subject to thermal vibrations.

Thus, its important to make sure that the crystal lattice does not drift as a whole during the simulation. Furthermore, it's important to make sure the crystal does not expand or contract, because that would also push atoms away from their original locations. To partially mitigate this type of problem, the WS modifier provides an option to rescale the expanded simulation cell, including all contained atomic coordinates, to the size of the reference cell.

Regarding your script:

The offset that you apply to the atomic positions has no effect, because it happens after the WS analysis is performed. You apply to modifiers, the WignerSeitzAnalysisModifier and the PythonScriptModifier, in that order. When the PythonScriptModifier moves the atoms, the WS modifier will not see it. You need to move the atoms with a modifier that is inserted before the WS modifier. If you don't want to use another PythonScriptModifier for this, you can also use a AffineTransformationModifier to translate all atoms by a uniform offset.
« Last Edit: September 26, 2017, 05:55:24 PM by Alexander Stukowski »

Alexander Stukowski

  • Administrator
  • Hero Member
  • *****
  • Posts: 638
Re: WS Analysis on system that is moving
« Reply #7 on: September 27, 2017, 01:26:57 PM »
Just as an addition to my previous description of the WS defect method: The OVITO user manual contains an illustration of the method, which might help understanding it:

http://ovito.org/manual/particles.modifiers.wigner_seitz_analysis.html