TST: Change r_sc test to look for vacancy issue (#265)

matscipy/dislocation.py

@@ -4513,7 +4513,7 @@ def _get_kink_quad_mask(self, ref_bulk, direction, precision=3):
 
 
         # 1st, 2nd, & 3rd nearest neighbour distances
-        # for each crytal structure
+        # for each crystal structure
         neigh_dists = {
         "bcc" : [0.88, 1.01, 1.42],
         "fcc" : [0.72, 1.01, 1.23],

tests/test_dislocation.py

@@ -21,6 +21,7 @@
 #
 import os
 import unittest
+from matscipy.neighbours import coordination
 import matscipytest
 import pytest
 import sys
@@ -30,7 +31,7 @@
 
 from ase.calculators.lammpslib import LAMMPSlib
 from matscipy.calculators.eam import EAM
-from matscipy.dislocation import DiamondGlide90degreePartial, get_elastic_constants
+from matscipy.dislocation import BCCEdge100Dislocation, DiamondGlide90degreePartial, get_elastic_constants
 from matscipy.calculators.manybody.explicit_forms.stillinger_weber import StillingerWeber,\
                                                                 Holland_Marder_PRL_80_746_Si
 from matscipy.calculators.manybody import Manybody
@@ -733,15 +734,59 @@ def test_cylinder_ovito_dxa(self, disloc, gen_bulk,
         # assert abs(err) < tol
 
     def test_displacement_r_sc(self, disloc, subtests):
+        '''
+        Test whether the r_sc parameter can fix the vacancy issue reported in https://github.com/libAtoms/matscipy/issues/265
+        
+        '''
+
+        allowed_disloc_classes = [
+            BCCEdge100Dislocation
+        ]
+
+        if disloc.__class__ not in allowed_disloc_classes:
+            self.skipTest()
+
         self.set_up_cls(disloc)      
 
         d = self.test_cls(self.alat, self.C11, self.C12, self.C44, symbol=self.symbol)
+
+        # This ref disloc should have a "vacancy" in it
+        # (Similar to #265)
+        ref_bulk, ref_disloc = d.build_cylinder(200.0, method=self.default_method, verbose=False)
+
+        # Get mask of all unfixed atoms (to exclude surfaces from coordination analysis)
+        outer_mask = ~ref_disloc.arrays["fix_mask"]
+
+        # Get mask of non-core atoms
+        p = ref_bulk.positions - ref_disloc.info["core_positions"][0]
+        inner_mask = np.linalg.norm(p[:, :2], axis=-1) > 15.0
+
+        full_mask = inner_mask * outer_mask
 
-        _, ref_disloc = d.build_cylinder(80.0, method=self.default_method, verbose=False)
+        # Overestimates of 1st nearest neighbour distances
+        # for each crystal structure (to allow for some elasticity)
+        neigh_dists = {
+        "bcc" : 0.89,
+        "fcc" : 0.73,
+        "diamond" : 0.45,
+        }
 
-        _, r_sc_disloc = d.build_cylinder(80.0, method=self.default_method, verbose=False, r_sc=60)
+        ref_coord, ref_counts = np.unique(coordination(ref_disloc, self.alat * neigh_dists[disloc.crystalstructure])[full_mask], return_counts=True)
 
-        np.testing.assert_almost_equal(ref_disloc.positions, r_sc_disloc.positions, decimal=2)
+        # This disloc should not
+        _, r_sc_disloc = d.build_cylinder(200.0, method=self.default_method, verbose=False, r_sc=100)
+
+        r_sc_coord, r_sc_counts = np.unique(coordination(r_sc_disloc, self.alat * neigh_dists[disloc.crystalstructure])[full_mask], return_counts=True)
+
+        bulk_coord = {
+            "bcc" : 8,
+            "fcc" : 12,
+            "diamond" : 4 
+        }
+
+        # Check that more atoms are bulk coordinated in r_sc_disloc than in ref_disloc
+        # (i.e. that the "vacancy" has been removed due to the more stable SC displacments)
+        assert ref_counts[ref_coord==bulk_coord[disloc.crystalstructure]] < r_sc_counts[r_sc_coord==bulk_coord[disloc.crystalstructure]]
 
     def test_glide_configs(self, disloc, subtests):  
         self.set_up_cls(disloc)