diff --git a/ase_sim.py b/ase_sim.py
index d270acef037011009d0e77b6976248e60407d56e..529e8b36fdd74d375ed59569bc9e84d8edc6b36d 100644
--- a/ase_sim.py
+++ b/ase_sim.py
@@ -3,6 +3,7 @@ import numpy as np
import math
import ase.lattice.cubic
+import ase.data
import ase.io
import ase.units
@@ -18,6 +19,7 @@ class AseSim:
self.atoms = None
self.dyn = None
self.run = None
+ self.tagged_elements = {}
def create_lattice(self, lattice_size, basis, empty_space, pot_filename, rand_ox=False):
self.atoms = ase.lattice.cubic.BodyCenteredCubic(basis, size=lattice_size, pbc=(True, True, False))
@@ -41,23 +43,22 @@ class AseSim:
nbor_list = NeighborList(nbor_cutoffs, skin=0.01, self_interaction=False, bothways=True)
nbor_list.update(self.atoms)
- pot_path = os.path.join('data', pot_filename)
-
- self.dyn = SajuhakDynamics(self.atoms, nbor_list, pot_path, rand_ox)
- self.run = self.dyn.run
- self.atoms.set_calculator(KokkoPotential(nbor_list, pot_path))
-
# initialize atom modifications
- symbols = self.dyn.atoms.get_chemical_symbols()
- positions = self.dyn.atoms.get_positions()
+ symbols = self.atoms.get_chemical_symbols()
+ positions = self.atoms.get_positions()
max_miller = lattice_size[2] * (1 - empty_space)
# replace atoms with empty spaces above a certain z position
for i, position in enumerate(positions):
- if position[2] > self.dyn.atoms.millerbasis[2][2]*max_miller:
+ if position[2] > self.atoms.millerbasis[2][2]*max_miller:
symbols[i] = 'X'
- self.dyn.atoms.set_chemical_symbols(symbols)
+ # store modifications and initialize dyn and calc
+ pot_path = os.path.join('data', pot_filename)
+ self.atoms.set_chemical_symbols(symbols)
+ self.dyn = SajuhakDynamics(self.atoms, nbor_list, pot_path, rand_ox)
+ self.run = self.dyn.run
+ self.atoms.set_calculator(KokkoPotential(nbor_list, pot_path))
# @profile # add -m memory_profiler to interpreter options
def dope(self, dopants, basis):
@@ -109,3 +110,26 @@ class AseSim:
if self.dyn.atoms.get_atomic_numbers()[nbor_idx] == 0:
empty_nbor_list[idx] += 1
self.dyn.atoms.set_array('empty_nbors', empty_nbor_list)
+
+ def tag_elements(self):
+ """
+ Change the tags of random elements in the lattice so they can be tracked
+
+ :return: None
+ """
+ tag_target = 5 # how many atoms to tag of each element
+ number_list = self.atoms.get_atomic_numbers()
+ element_list = np.unique(number_list)
+ tag_list = self.atoms.get_tags()
+ tagged_counts = {} # dict of amount of tagged atoms per element
+ tag_index = 1
+ while sum(tagged_counts.values()) < tag_target * len(element_list):
+ rand_idx = np.random.randint(len(self.atoms))
+ rand_number = number_list[rand_idx]
+ tagged_count = tagged_counts.get(rand_number, 0)
+ if not tagged_count >= tag_target:
+ tag_list[rand_idx] = tag_index
+ self.tagged_elements[tag_index] = rand_number
+ tag_index += 1
+ tagged_counts[rand_number] = tagged_count + 1
+ self.atoms.set_tags(tag_list)
diff --git a/sajuhak_dyn.py b/sajuhak_dyn.py
index 935d16d05898c3d46165f351249b3f18bc67791c..4838c5ad6a9790dd5e8c5072ed44764d1493ad64 100644
--- a/sajuhak_dyn.py
+++ b/sajuhak_dyn.py
@@ -97,7 +97,7 @@ class SajuhakDynamics(Dynamics):
self.nbor_list = nbor_list
# self.nbor_list.update(self.atoms)
self.rand_ox = rand_ox
- self.diff_const_dict = {}
+ self.tag_count_list = {}
def run(self, kT, n_steps=100):
for i in range(n_steps):
@@ -111,38 +111,26 @@ class SajuhakDynamics(Dynamics):
rand_idx = np.random.randint(len(self.atoms))
rand_number = self.atoms.get_atomic_numbers()[rand_idx]
rand_symbol = ase.atoms.chemical_symbols[rand_number]
- if rand_symbol not in self.diff_const_dict.keys():
- self.diff_const_dict[rand_symbol] = [0, 0]
- self.diff_const_dict[rand_symbol][1] += 1
nbor_indices = self.nbor_list.neighbors[rand_idx]
nbor_idx = np.random.choice(nbor_indices)
nbor_number = self.atoms.get_atomic_numbers()[nbor_idx]
nbor_symbol = ase.atoms.chemical_symbols[nbor_number]
- if nbor_symbol not in self.diff_const_dict.keys():
- self.diff_const_dict[nbor_symbol] = [0, 0]
- self.diff_const_dict[nbor_symbol][1] += 1
if rand_number == nbor_number:
- # if the two atoms are the same, nothing changes, remove from total exchanges
- self.diff_const_dict[rand_symbol][1] -= 1
- self.diff_const_dict[nbor_symbol][1] -= 1
+ # if the two atoms are the same, nothing changes
return
# if either point is empty, initiate random (de)oxidation
elif self.rand_ox and (rand_number == 0 or nbor_number == 0):
# if either point is not oxygen, initiate random oxidation
if rand_number != 8 and nbor_number != 8:
if self.random_oxidation(kT, rand_idx, nbor_idx):
- # if the oxidation was successful, don't swap these atoms yet, remove from total exchanges
- self.diff_const_dict[rand_symbol][1] -= 1
- self.diff_const_dict[nbor_symbol][1] -= 1
+ # if the oxidation was successful, don't swap these atoms yet
return
# else initiate random deoxidation
else:
if self.random_deoxidation(kT, rand_idx, nbor_idx):
- # if the deoxidation was successful, don't swap these atoms yet, remove from total exchanges
- self.diff_const_dict[rand_symbol][1] -= 1
- self.diff_const_dict[nbor_symbol][1] -= 1
+ # if the deoxidation was successful, don't swap these atoms yet
return
energy_start = self.calc_energy_at_points([rand_idx, nbor_idx])
@@ -154,8 +142,14 @@ class SajuhakDynamics(Dynamics):
if not np.random.random() < exp(-delta_e / kT):
self.swap_atoms(nbor_idx, rand_idx)
else:
- self.diff_const_dict[rand_symbol][0] += 1
- self.diff_const_dict[nbor_symbol][0] += 1
+ # swap tags
+ tag_list = self.atoms.get_tags()
+ rand_tag = tag_list[rand_idx]
+ self.tag_count_list[rand_tag] = self.tag_count_list.get(rand_tag, 0) + 1
+ nbor_tag = tag_list[nbor_idx]
+ self.tag_count_list[nbor_tag] = self.tag_count_list.get(nbor_tag, 0) + 1
+ tag_list[rand_idx], tag_list[nbor_idx] = tag_list[nbor_idx], tag_list[rand_idx]
+ self.atoms.set_tags(tag_list)
def swap_atoms(self, index1, index2):
atom1_number = self.atoms[index1].number
diff --git a/sim_controller.py b/sim_controller.py
index a146eed8e85343267a7bd22c661634d37d39a5b3..7645ba4a8e79195942534652958b365b9e2b956b 100644
--- a/sim_controller.py
+++ b/sim_controller.py
@@ -13,6 +13,7 @@ from collections import OrderedDict
# import cProfile
import ase_sim
+from ase.data import chemical_symbols
import vis
@@ -175,6 +176,7 @@ class SimController(object):
self.sim.create_lattice(**self.lattice_params, rand_ox=self.sim_params['rand_ox'])
self.bulk_counts = self.sim.dope(self.bulk_dopants, self.sim.basis)
self.empty_space_counts = self.sim.dope(self.empty_space_dopants, 'X')
+ self.sim.tag_elements()
def run(self):
if type(self.sim_params['kT']) is list:
@@ -223,23 +225,22 @@ class SimController(object):
lattice_counts = bulk_counts_header + bulk_counts_str + vacuum_counts_header + vacuum_counts_str
- diff_const_header = "Exchanges per element (successful/total (percentage)):\n"
- diff_const_line = " {symbol}: {successful}/{total} ({perc:.2f}%)\n"
+ diff_const_header = "Distances traveled by tagged atoms (average, list):\n"
+ diff_const_line = " {symbol}: {average}, {list}\n"
diff_const_str = ""
- # tot_exchanges = self.sim_params['n_exchanges']
- # if type(tot_exchanges) is list:
- # tot_exchanges = sum(tot_exchanges)
-
- diff_const_dict_items_sorted_by_total = sorted(self.sim.dyn.diff_const_dict.items(),
- key=lambda item: item[1][1])
- for symbol, success_tot in diff_const_dict_items_sorted_by_total:
- if symbol == 'X':
- symbol = 'Vac'
- diff_const_str += diff_const_line.format(symbol=symbol,
- successful=success_tot[0],
- total=success_tot[1],
- perc=success_tot[0] / success_tot[1] * 100)
+ counts_by_element = {}
+ for index, number in self.sim.tagged_elements.items():
+ count = self.sim.dyn.tag_count_list.get(index, 0)
+ if number not in counts_by_element.keys():
+ counts_by_element[number] = []
+ counts_by_element[number].append(count)
+
+ for number, counts in sorted(counts_by_element.items()):
+ if number != 0:
+ diff_const_str += diff_const_line.format(symbol=chemical_symbols[number],
+ average=np.mean(counts),
+ list=counts)
diff_const = diff_const_header + diff_const_str
diff --git a/vis.py b/vis.py
index 73827315bfcb33637f611639a5478815da4ad10c..bf9f3f73d70c11ec7eabd213cf2f2b534552b45f 100644
--- a/vis.py
+++ b/vis.py
@@ -26,15 +26,17 @@ def read_xyz_frame(xyz_file_in):
list_properties = []
comment = xyz_file_in.readline()
info = ase.io.extxyz.key_val_str_to_dict(comment)
- # properties, names, dtype, convs = ase.io.extxyz.parse_properties(info['Properties'])
+ _, names, _, _ = ase.io.extxyz.parse_properties(info['Properties'])
# read the next n_sites lines, the index is not used
for _ in range(n_sites):
- # an extxyz line is "element x y z properties"
line = xyz_file_in.readline()
- if info['Properties'].find('empty_nbors:I:1') > info['Properties'].find('Z:I:1'):
- element, x, y, z, _, empty_nbors = line.split()
- else:
- element, x, y, z, empty_nbors, _ = line.split()
+ values = line.split()
+ element, x, y, z = values[:4]
+ values = values[4:]
+ row = {}
+ for i, name in enumerate(names[2:]):
+ row[name] = values[i]
+ empty_nbors = row['empty_nbors']
# add to the list of particles
list_properties.append((element, float(x), float(y), float(z), int(empty_nbors)))
# return the list of particles