Purpose

This script facilitates the calculation of the velocity of a dislocation or other moving defect. The script requires only a LAMMPS dump file containing atom positions and will output a data file of timestep and defect position.

The script makes several assumptions.

• The bulk structure is the one shared by the most atoms.
• The defect is in the center of the Z direction of the cell and can be identified by not sharing the bulk structure.
• The defect is moving in the X direction.
• The cell is periodic in the X direction.

The script also requires version 2.7 or newer of Ovito.

Usage

There are two ways to use the script, within the graphical Ovito program, or from the command line using the Ovitos scripting utility.

Using Ovitos Scripting Utility (Recommended)

From the command line, simply run the command below. You may need to specify explicitly the path to the ovitos executable in your Ovito distribution.

ovitos dis_mobility.py dump.atom.positions

Using the Ovito GUI

In the gui, you will need to add four modifiers such that they are in the following order (top to bottom) in the modifier list.

• Python Script
• Delete Selected Particles
• Expression Select
• Polyhedral Template Matching

Polyhedral Template Matching

The Polyhedral Template Matching modifier will determine the crystal structure on a per atom basis. It requires no setup for this case.

Expression Select

We will use this modifier to select all the atoms that are not a part of the defect that we are concerned about. In the expression box, add the conditional statement

StructureType == BULK || Position.Y < BELOW_MID || Position.Y > ABOVE_MID

Where BULK is 1 if your material has an FCC bulk structure, 2 for HCP, or 3 for BCC.BELOW_MID and ABOVE_MID should be a distance above and below the defect, for example, 30 Angstroms below and above the center of the cell in the Z direction, respectively.

Delete Selected Particles

This modifier will delete all selected particles, as the name says. It requires no setup.

Python Script

This modifier will load the code from the script to create the defect position data file. You will need to make sure that the source file below is somewhere on your PYTHONPATH. Then, add the following to the script window, keeping the default import lines, but replacing the default modify function definition.

from dis_mobility import get_dislocation_position
modify = get_dislocation_position

Source

# Calculate the speed of a dislocation
# Bradley Huddleston, February 15, 2017
#

# To use this within Ovito in the Python Script modifier,
#   use the following lines in triple quotes as the python script.
#   Replace SIMPATH with your simulation directory, and '99' with 
#   the number of timesteps in your output file.
"""
from ovito.data import *
import os
os.chdir('SIMPATH')
import dis_mobility as dm

modify = dm.get_dislocation_position
dm.initialize_data(99)
"""


from ovito import *
import numpy as np
import os

assert(version[0] >= 2 and version[1] >= 7)


def get_dislocation_position(frame, input, output):

	x_extent = output.cell.matrix[0,0]
	positions = output.particle_properties.position.marray #.copy()
	
	# Eliminate all periodic images
	positions[(positions[:,0] > x_extent),0] = positions[(positions[:,0] > x_extent),0] - x_extent
	# Edit positions to put dislocation all on same side 
	if frame == 0:
		global core_spread
		core_spread = np.max(positions[:,0]) - np.min(positions[:,0])
		if core_spread > 0.5*x_extent:
			core_spread = None 
	elif core_spread:
		ref = positions[0,0]
		for i,p in enumerate(positions):
			if i == 0:
				continue 
			if (p[0] - ref) > core_spread*1.3:
				p[0] -= x_extent
			ref = np.mean(positions[0:(i+1),0])
	dislocation = np.mean(positions,0)[0]
	timestep = output.attributes['Timestep']
	
	while (pos[frame - 1] - dislocation) > 0.25*x_extent:
		dislocation += x_extent

	with open('posVStime.txt','a+') as f:
		f.write('{}\t{}\n'.format(timestep, dislocation))
	
	time[frame] = timestep
	pos[frame] = dislocation	


def initialize_data(num):
	global time, pos, core_spread
	time = np.zeros(num)
	pos = np.zeros(num)
	core_spread = None
	
	if os.path.exists('posVStime.txt'):
		data = np.loadtxt('posVStime.txt')
		time.put(range(len(data[:,0])),data[:,0])
		pos.put(range(len(data[:,1])),data[:,1])
		

# If run from the command line setup required modifiers
if __name__ == "__main__":

	from sys import argv

	if len(argv) < 2:
		raise	Exception("Filename required as an argument.")

	filename = argv[-1]
	
	# Parse input arguments
	# Default values
	if len(argv) > 2:
		keyargs = argv[1:-1]
		i = 0
		while (i < len(keyargs)):
			key = keyargs[i].split('-')[1]
			i += 1
			if (key is "h") or (key is "help"):
				print("Usage: ovitos dis_mobility.py [OPTIONS] [FILENAME]\n\t-h,help\tThis help message\n")
			else:
				raise Exception("Unknown option: " + key)

	# ### Setup ovitos ###
	# Import file
	if len([s for s in filename if s == '*']) > 0:
		node = io.import_file(filename) 
	else:
		node = io.import_file(filename, multiple_frames=True)

	# Create CNA or PTM modifier
	#	if version[1] = 7:
		#structure = modifiers.CommonNeighborAnalysisModifier()
		#elif version[1] > 7:
	structure = modifiers.PolyhedralTemplateMatchingModifier()
	node.modifiers.append(structure)
	node.compute()
	
	# Estimate bulk structure
	nfcc = node.output.attributes['PolyhedralTemplateMatching.counts.FCC']
	nbcc = node.output.attributes['PolyhedralTemplateMatching.counts.BCC']
	nhcp = node.output.attributes['PolyhedralTemplateMatching.counts.HCP']
	struct_ID = np.argmax([nfcc,nhcp,nbcc]) + 1
	
	# Add expression select
	y_extent = node.output.cell.matrix[1,1]
	bulk = modifiers.SelectExpressionModifier(
		expression = 'StructureType == {} || Position.Y < {} || Position.Y > {}'.format(
			struct_ID,
			y_extent/2.0 - 30.0,
			y_extent/2.0 + 30.0))
	node.modifiers.append(bulk)
	
	# Add delete selected
	delete = modifiers.DeleteSelectedParticlesModifier()
	node.modifiers.append(delete)
	
	numframes = node.source.num_frames
	frames = range(numframes)

	initialize_data(numframes)

	for frame in frames:
		if frame < numframes:
			dataset.anim.current_frame = frame # Advance frame. 
		else:
			continue
	
		node.compute()
		
		get_dislocation_position(frame, node.output, node.output)