Skip to content

Commit

Permalink
TST: Gamma Surface tests
Browse files Browse the repository at this point in the history
  • Loading branch information
thomas-rocke committed Oct 5, 2023
1 parent 1349682 commit 4977379
Show file tree
Hide file tree
Showing 2 changed files with 33 additions and 126 deletions.
16 changes: 8 additions & 8 deletions matscipy/gamma_surface.py
Original file line number Diff line number Diff line change
Expand Up @@ -118,7 +118,7 @@ def _y_dir_search(self, d1):
# No nice integer vector found!
raise RuntimeError(f"Could not automatically find an integer basis from basis vector {d1}")

def generate_images(self, nx, ny, z_replications=1, atom_offset=0.0, cell_strain=0.0, vacuum=0.0, vacuum_offset=0.0, path_limsx=[0.0, 1.0], path_limsy=None):
def generate_images(self, nx, ny, z_replications=1, atom_offset=0.0, cell_strain=0.0, vacuum=0.0, vacuum_offset=0.0, path_xlims=[0.0, 1.0], path_ylims=None):
'''
Generate gamma surface images on an (nx, ny) grid
Expand All @@ -138,19 +138,19 @@ def generate_images(self, nx, ny, z_replications=1, atom_offset=0.0, cell_strain
Offset (in A) applied to the position of the vacuum layer in the cell
The position of the vacuum layer is given by:
vac_pos = self.cut_at.cell[2, 2] * (1 + cell_strain) / 2 + vacuum_offset
path_limsx: list|array of floats
path_xlims: list|array of floats
Limits (in fractional coordinates) of the stacking fault path in the x direction
path_limsy: list|array of floats
path_ylims: list|array of floats
Limits (in fractional coordinates) of the stacking fault path in the x direction
If not supplied, will be set to path_limsx
If not supplied, will be set to path_xlims
'''

if path_limsy is None:
y_lims = path_limsx
if path_ylims is None:
y_lims = path_xlims
else:
y_lims = path_limsy
y_lims = path_ylims

x_lims = path_limsx
x_lims = path_xlims

self.nx = nx
self.ny = ny
Expand Down
143 changes: 25 additions & 118 deletions tests/test_gamma_surface.py
Original file line number Diff line number Diff line change
Expand Up @@ -9,127 +9,34 @@

class GammaSurfaceTest(matscipytest.MatSciPyTestCase):
def setUp(self) -> None:
self.fmax = 1e-4
self.at0 = Diamond('Si', latticeconstant=5.43)

# Si testing framework SW model string
# https://github.com/libAtoms/silicon-testing-framework/blob/master/models/SW/model.py
self.model = Potential('IP SW', param_str="""
<SW_params n_types="1">
<per_type_data type="1" atomic_num="14" />
<per_pair_data atnum_i="14" atnum_j="14" AA="7.049556277" BB="0.6022245584"
p="4" q="0" a="1.80" sigma="2.0951" eps="2.1675" />
<per_triplet_data atnum_c="14" atnum_j="14" atnum_k="14"
lambda="21.0" gamma="1.20" eps="2.1675" />
</SW_params>
""")
self.fmax = 1e-3
self.at0 = Diamond('Si', latticeconstant=5.43)

# Si testing framework SW model string
# https://github.com/libAtoms/silicon-testing-framework/blob/master/models/SW/model.py
self.model = Potential('IP SW', param_str="""
<SW_params n_types="1">
<per_type_data type="1" atomic_num="14" />
<per_pair_data atnum_i="14" atnum_j="14" AA="7.049556277" BB="0.6022245584"
p="4" q="0" a="1.80" sigma="2.0951" eps="2.1675" />
<per_triplet_data atnum_c="14" atnum_j="14" atnum_k="14"
lambda="21.0" gamma="1.20" eps="2.1675" />
</SW_params>
""")

def test_gamma_surface(self):
# Test functionality of gamma surface code + using model to relax

surface = GammaSurface(self.at0, np.array([0, 0, 1], dtype=int), np.array([1, 0, 0], dtype=int))
surface.generate_images(7, 7)
surface.relax_images(self.model, self.fmax)
Es = surface.get_surface_energies(self.model)
surface = GammaSurface(self.at0, np.array([0, 0, 1], dtype=int), np.array([1, 0, 0], dtype=int))
surface.generate_images(3, 3)
surface.relax_images(self.model, self.fmax)
Es = surface.get_surface_energies(self.model)

assert np.min(Es) >= 0.0
print(np.max(Es))

assert np.min(Es) >= 0.0
assert np.allclose(np.max(Es), 0.1946074084327025)

def test_stacking_fault(self):
import matplotlib.pyplot as plt
from ase.io import write
# Test that stacking fault reproduces Si results from Si testing framework for SW potential
surface = StackingFault(self.at0, np.array([0, 0, 1]), np.array([1, 1, 0]), compress=True)
surface.generate_images(21, z_replications=1)
surface.relax_images(self.model, self.fmax)
Es = surface.get_surface_energies(self.model)

surface.plot_gamma_surface()
plt.savefig("Stack_Compressed.png")

Es1 = Es.copy()[0, :]
l1 = len(surface.images[0])

surf1 = surface.surface_area

end_E = surface.images[4].get_potential_energy()
# write("Surface_structs_compressed.xyz", surface.images)

xd1 = surface.x_disp
yd1 = surface.y_disp
sep1 = surface.surface_separation

cell1 = surface.images[0].cell
surface1 = surface

surface = StackingFault(self.at0, np.array([0, 0, 1]), np.array([1, 1, 0]), compress=False)
surface.generate_images(21, z_replications=1)
surface.relax_images(self.model, self.fmax)
Es = surface.get_surface_energies(self.model)[0, :]
l2 = len(surface.images[0])

surface.plot_gamma_surface()
plt.savefig("Stack_Uncompressed.png")




surface = GammaSurface(self.at0, np.array([0, 0, 1]), np.array([1, 1, 0]), compress=True)
surface.generate_images(21, 21, z_replications=1)
surface.relax_images(self.model, self.fmax)
Es = surface.get_surface_energies(self.model)[0, :]
l2 = len(surface.images[0])

surface.plot_gamma_surface()
plt.savefig("Gamma_Compressed.png")

surface = GammaSurface(self.at0, np.array([0, 0, 1]), np.array([1, 1, 0]), compress=False)
surface.generate_images(21, 21, z_replications=1)
surface.relax_images(self.model, self.fmax)
Es = surface.get_surface_energies(self.model)[0, :]
l2 = len(surface.images[0])

surface.plot_gamma_surface()
plt.savefig("Gamma_Uncompressed.png")


print(xd1)
print(yd1)

print()

print(surface.x_disp)
print(surface.y_disp)

print()
print(surface1.mapping)

# print(Es1.shape, Es.shape)

# print(Es)
# print(Es1)

# print(l1, l2)
# print(cell1[:, :])
# print(surface.images[0].cell[:, :])

# write("Surface_structs_uncompressed.xyz", surface.images)

print(np.max(Es1), np.max(Es))

print(np.max(Es1) * surf1, np.max(Es) * surface.surface_area)


tst = GammaSurfaceTest()
tst.setUp()
#tst.test_gamma_surface()
tst.test_stacking_fault()

## Si paper gives 2.17 J/m**2 == 0.1354 eV/A**2 for (0 0 1)
## compressed result = 0.16738092
## uncompressed = 0.07231819518
surface = StackingFault(self.at0, np.array([1, 1, 0]), np.array([-1, 1, 0]))
surface.generate_images(9, z_replications=1, path_ylims=[0, 0.5])
surface.relax_images(self.model, self.fmax)
Es = surface.get_surface_energies(self.model)[0, :]

#Smaller
## Compressed = 0.1673890
## Uncompressed = 0.723372
assert np.allclose(np.max(Es), 0.16738903450149265)

0 comments on commit 4977379

Please sign in to comment.