tag atoms and count their exchanges

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)

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

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
 

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