Add diffusion constant testing

sajuhak_dyn.py

@@ -97,6 +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 = {}
 
     def run(self, kT, n_steps=100):
         for i in range(n_steps):
@@ -109,9 +110,18 @@ class SajuhakDynamics(Dynamics):
         # choose a random point and one of its neighbors
         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
 
         # check what we chose
         if rand_number == nbor_number:
@@ -125,19 +135,18 @@ class SajuhakDynamics(Dynamics):
             else:
                 if self.random_deoxidation(kT, rand_idx, nbor_idx):
                     return
-        # start_energy = self.atoms.get_potential_energy()
-        start_energy = self.calc_energy_at_points([rand_idx, nbor_idx])
 
+        start_energy = self.calc_energy_at_points([rand_idx, nbor_idx])
         self.swap_atoms(rand_idx, nbor_idx)
-        # self.nbor_list.update(self.atoms)
-
-        # end_energy = self.atoms.get_potential_energy()
         end_energy = self.calc_energy_at_points([rand_idx, nbor_idx])
+
         delta_e = end_energy - start_energy
         # if the energy of the system increased, revert the change
         if not np.random.random() < exp(-delta_e / kT):
             self.swap_atoms(nbor_idx, rand_idx)
-            # self.nbor_list.update(self.atoms)
+        else:
+            self.diff_const_dict[rand_symbol][0] += 1
+            self.diff_const_dict[nbor_symbol][0] += 1
 
     def swap_atoms(self, index1, index2):
         atom1_number = self.atoms[index1].number