"""Modules for making crystallographic plane surfaces."""
from jarvis.core.atoms import Atoms
from jarvis.core.utils import ext_gcd
import numpy as np
from jarvis.analysis.structure.spacegroup import Spacegroup3D
from numpy.linalg import norm
from numpy import gcd
from collections import OrderedDict
[docs]def wulff_normals(miller_indices=[], surface_energies=[]):
"""Obtain Wulff Normals.
Args:
miller_indices : Miller indices
surface_energies : corresponding surface energies
Returns: Surface normals
"""
max_s = min(surface_energies)
normals = []
for i, j in zip(miller_indices, surface_energies):
normal = j * np.linalg.norm(i) / float(max_s)
normals.append([normal, i])
return normals
[docs]class Surface(object):
"""Get surface object of arbitrary atoms object and miller index."""
def __init__(
self,
atoms=None,
indices=[0, 0, 1],
layers=3,
thickness=25,
vacuum=18.0,
tol=1e-10,
from_conventional_structure=True,
):
"""Initialize the class.
Args:
atoms: jarvis.core.Atoms object
indices: Miller indices
layers: Number of surface layers
thickness: Provide thickness instead of layers
vacuum: vacuum padding
tol: tolerance during dot product
from_conventional_structure: whether to use the conv. atoms
"""
self.indices = np.array(indices)
self.from_conventional_structure = from_conventional_structure
if self.from_conventional_structure:
self.atoms = Spacegroup3D(atoms).conventional_standard_structure
else:
self.atoms = atoms
self.tol = tol
self.vacuum = vacuum
self.layers = layers
self.thickness = thickness
# Note thickness overwrites layers
[docs] def to_dict(self):
"""Convert to a dictionary."""
d = OrderedDict()
d["atoms"] = self.atoms.to_dict()
d["indices"] = self.indices
d["tol"] = self.tol
d["vacuum"] = self.vacuum
d["layers"] = self.layers
d["from_conventional_structure"] = self.from_conventional_structure
return d
[docs] @classmethod
def from_dict(self, d={}):
"""Construct class from a dictionary."""
return Surface(
atoms=Atoms.from_dict(d["atoms"]),
indices=d["indices"],
tol=d["tol"],
vacuum=d["vacuum"],
layers=d["layers"],
from_conventional_structure=d["from_conventional_structure"],
)
[docs] def make_surface(self):
"""Generate specified surface. Modified from ase package."""
atoms = self.atoms
h_index, k_index, l_index = self.indices
h0, k0, l0 = self.indices == 0
if h0 and k0 or h0 and l0 or k0 and l0: # if two indices are zero
if not h0:
c1, c2, c3 = [(0, 1, 0), (0, 0, 1), (1, 0, 0)]
if not k0:
c1, c2, c3 = [(0, 0, 1), (1, 0, 0), (0, 1, 0)]
if not l0:
c1, c2, c3 = [(1, 0, 0), (0, 1, 0), (0, 0, 1)]
else:
p, q = ext_gcd(k_index, l_index)
a1, a2, a3 = self.atoms.lattice_mat # .lat_lengths()
# constants describing the dot product of basis c1 and c2:
# dot(c1,c2) = k1+i*k2, i in Z
k1 = np.dot(
p * (k_index * a1 - h_index * a2)
+ q * (l_index * a1 - h_index * a3),
l_index * a2 - k_index * a3,
)
k2 = np.dot(
l_index * (k_index * a1 - h_index * a2)
- k_index * (l_index * a1 - h_index * a3),
l_index * a2 - k_index * a3,
)
if abs(k2) > self.tol:
i = -int(round(k1 / k2))
p, q = p + i * l_index, q - i * k_index
a, b = ext_gcd(p * k_index + q * l_index, h_index)
c1 = (p * k_index + q * l_index, -p * h_index, -q * h_index)
c2 = np.array((0, l_index, -k_index)) // abs(gcd(l_index, k_index))
c3 = (b, a * p, a * q)
lattice = atoms.lattice_mat # .lat_lengths()
basis = np.array([c1, c2, c3])
scaled = np.linalg.solve(basis.T, np.array(atoms.frac_coords).T).T
scaled -= np.floor(scaled + self.tol)
new_coords = scaled
tmp_cell = np.dot(basis, lattice)
M = np.linalg.solve(lattice, tmp_cell)
cart_coords = np.dot(scaled, lattice)
new_coords = np.dot(cart_coords, M)
new_atoms = Atoms(
lattice_mat=tmp_cell,
coords=new_coords,
elements=atoms.elements,
cartesian=True,
)
if self.thickness is not None and (self.thickness) > 0:
self.layers = int(self.thickness / new_atoms.lattice.c) + 1
# dims=get_supercell_dims(new_atoms,enforce_c_size=self.thickness)
# print ('dims=',dims,self.layers)
# surf_atoms = new_atoms.make_supercell_matrix([1, 1, dims[2]])
# print('self.layers',self.layers,self.thickness,new_atoms.lattice.c)
surf_atoms = new_atoms.make_supercell_matrix([1, 1, self.layers])
# print("supercell_cart_coords", surf_atoms.frac_coords)
new_lat = surf_atoms.lattice_mat # lat_lengths()
a1 = new_lat[0]
a2 = new_lat[1]
a3 = new_lat[2]
new_lat = np.array(
[
a1,
a2,
np.cross(a1, a2)
* np.dot(a3, np.cross(a1, a2))
/ norm(np.cross(a1, a2)) ** 2,
]
)
a1 = new_lat[0]
a2 = new_lat[1]
a3 = new_lat[2]
# print("a1,a2,a3", new_lat)
latest_lat = np.array(
[
(np.linalg.norm(a1), 0, 0),
(
np.dot(a1, a2) / np.linalg.norm(a1),
np.sqrt(
np.linalg.norm(a2) ** 2
- (np.dot(a1, a2) / np.linalg.norm(a1)) ** 2
),
0,
),
(0, 0, np.linalg.norm(a3)),
]
)
M = np.linalg.solve(new_lat, latest_lat)
new_cart_coords = surf_atoms.cart_coords # np.dot(scaled,lattice)
new_coords = np.dot(new_cart_coords, M)
new_atoms = Atoms(
lattice_mat=latest_lat,
elements=surf_atoms.elements,
coords=new_coords,
cartesian=True,
).center_around_origin()
frac_coords = new_atoms.frac_coords
frac_coords[:] = frac_coords[:] % 1
new_atoms = Atoms(
lattice_mat=latest_lat,
elements=surf_atoms.elements,
coords=frac_coords,
cartesian=False,
)
new_lat = new_atoms.lattice_mat
new_cart_coords = new_atoms.cart_coords
elements = new_atoms.elements
new_lat[2][2] = new_lat[2][2] + self.vacuum
with_vacuum_atoms = Atoms(
lattice_mat=new_lat,
elements=elements,
coords=new_cart_coords,
cartesian=True,
)
# new_atoms.center()
# print (with_vacuum_atoms)
return with_vacuum_atoms
"""
if __name__ == "__main__":
box = [[2.715, 2.715, 0], [0, 2.715, 2.715], [2.715, 0, 2.715]]
coords = [[0, 0, 0], [0.25, 0.25, 0.25]]
elements = ["Si", "Si"]
Si = Atoms(lattice_mat=box, coords=coords, elements=elements)
Surface(atoms=Si, indices=[1, 1, 1]).make_surface()
su = [
0.8582640971273426,
0.9334963319196496,
0.9360461382184894,
0.9419095687284446,
0.9802042233627004,
0.9875446840480956,
1.0120634294466684,
1.0126231880823566,
1.0241538763302507,
1.0315901848682645,
1.0318271257831195,
1.0331286888257398,
1.0344297141291043,
1.0388709097092674,
1.040277640596931,
1.042494119906149,
1.04453679643896,
1.0450598648770613,
1.045076130339553,
1.0469310544190567,
1.0491015867538047,
1.0495494553198788,
1.0534717916897114,
1.0535201391639715,
1.054233162444997,
1.0579157863887743,
1.0595676718662346,
1.0601381085497692,
1.109580394178689,
]
ml = [
[0, 0, 1],
[2, 0, 3],
[2, 0, 1],
[1, 0, 1],
[3, 0, 2],
[1, 0, 3],
[3, 1, 1],
[3, 0, 1],
[3, 1, 3],
[3, -1, 1],
[3, 1, 0],
[3, 2, 1],
[3, 3, 1],
[1, 0, 0],
[2, 2, 1],
[3, -1, 3],
[3, -1, 2],
[3, 3, 2],
[3, 2, 2],
[2, -1, 3],
[3, 2, 0],
[3, 2, 3],
[1, 1, 1],
[1, 0, 2],
[3, 1, 2],
[2, -1, 2],
[3, -1, 0],
[2, 2, 3],
[1, 1, 0],
]
nm = wulff_normals(miller_indices=ml, surface_energies=su)
print(nm)
from jarvis.core.lattice import Lattice
lat = Lattice([[4.05, 0, 0], [0, 4.05, 0], [0, 0, 4.05]])
pmg_wulff = WulffShape(lat, ml, su)
print(pmg_wulff.facets)
"""