# Introduction
In this note, I will discuss a semi-automated process of constructing nanowires grown along arbitrary directions.
The only manual part of the process is figuring out the transformation matrix so that the periodic direction (usually $z$) is oriented along the growth direction (if you know how to automate this, please let me know!). Everything from that point forward is automated using `pymatgen`.
I will be using the antiperovskite Sr$_3$PbO as an example. Note that this system is cubic, so some of the assumptions made here do not hold for non-cubic systems.
> [!TIP]
> If your system is hexagonal, you can find a transformation matrix to make it orthorhombic.
>
> However, for orthorhombic cells, it is not quite as straightforward to determinate the transformation matrix if you have some weird growth direction like $[\bar{2}11]$ that I use below.
# The Transformation Matrix
In this section, we will compute the transformation matrix that allows us to construct a "primitive" nanowire cell whose $c$ axis pointing along the $[\bar{2}11]$ direction of the primitive cell.
## More Atoms? More Atoms.
We start by adding some more unit cells using the `boundary` function, do not use the transform function for this specific case. I started with the boundaries $[-5,4]\times[-5,4]\times[-5,4]$. Make sure you delete any bonds.
## Plane 1
We want the nanowire supercell to be periodic in the $[\bar{2}11]$ direction, so we can place a $(\bar{2}11)$ plane using `edit -> lattice planes`
It doesn't matter where you put this plane for now, I just picked $0 \times d$ . At this point, you'd normally rotate the structure until the normal to the new plane points up, and rotate along that axis until you have some nicely aligned atoms. For example, see the picture below.
> [!WARNING]
> For non-cubic cells, you generally won't find such a neat row of atoms perpendicular to the plane. You'll have to get a little creative, and the resultant nanowire cell might not be orthorhombic.
![[Sr3PbO_primnw_planes 5.png]]
## Plane 2
The next step is to figure out the miller indices of planes where we'll trim the system laterally. The highlighted atoms in the figure above (including the ones behind them) show an example of a plane orthogonal to the $(\bar{2}11)$ plane (technically, the planes need not be orthogonal, but this allows for the construction of an orthorhombic cell).
### Getting the miller indices automatically
With the atoms selected, add a plane and use the "calculate the best plane for the selected atoms". This will give you the miller indices right away, in this case, $(01\bar{1})$ (you might get something that slightly deviates from these numbers, just round it).
### Getting the miller indices manually
If you want to do this by hand for whatever reason, pick the coordinates of any three non-collinear atoms on the plane:
$
\begin{align}
\text{Pb}_1&: (-3, -2, 2) \equiv P, \\
\text{Pb}_2&: (1, -3, 1) \equiv Q, \\
\text{Pb}_3&: (-1, -3, 1) \equiv R.
\end{align}
$
Now we have the coordinates of 3 points, and we can easily obtain the normal vector to the plane as $\vec{v} = \overrightarrow{PQ} \times \overrightarrow{PR}$ .
Here's some code that does that, because I'm lazy.
```python
import numpy as np
P = np.array([2.5, 0.0, 6.5])
Q = np.array([3.0, 0.0, 6.0])
R = np.array([3.0, 0.5, 5.5])
PQ = Q - P
PR = R - P
v = np.cross(PQ, PR)
```
We end up with:
$
\mathbf{v} = 2\times(0,1,-1)
$
Getting rid of the overall factor, this means the plane we're seeking has a normal vector $[01\bar{1}]$. Since the system is cubic, the plane itself is also $(01\bar{1})$.
> [!WARNING]
> If your system is not cubic, the normal to the plane $(hkl)$ is not given by $[hkl]$. It's not awfully hard to figure it out the miller indices, however.
## Plane 3
To get the normal vector to the final constraining plane (and thus the miller indices of the plane itself, since the system is cubic), we can simply take the cross product of the normal vectors of the two other planes, and we get $(-2, -2, -2)$ . This means the miller indices of the last plane are $(111)$.
> [!TIP]
> Now that we know the miller indices of all three constraining planes, it's easy to use the orientation dialogue to get the perfect view.
>
> First select "Project along normal to $(hkl)
quot;. The upward vector is easy, we want the new $z$ direction to be orthogonal to the $(\bar{2}11)$ plane, and thus parallel to $[\bar{2}11]$, so we enter $(u,v,w) = (-2,1,1)$.
>
> The projection vector is a little tricky. Essentially, it is the vector that points out of the screen, orthogonal to the upward vector. We already know two good options, $[111]$ and $[01\bar{1}]$.
## Constructing the Transformation Matrix
### The Easy Way
Since the system is cubic, the transformation matrix is simply:
$
T = \begin{pmatrix}
1 & 1 & 1 \\
0 & -1 & 1 \\
2 & -1 & -1
\end{pmatrix}
$
The first and second rows can be switched, the two resultant cells are simply related by a rotation. You can also multiply any row by $-1$ and the result is still the same, again related by a rotation.
Note that $\det{T} = 6$, i.e., the new cell consists of 6 primitive unit cells.
### The Hard Way
This method is quite involved and (naturally) gives the same answer. I'm including it for informational purposes.
#### Adding Planes
Let's start by adding 2 of each of the 3 planes we have using the `lattice planes` dialogue. Now we need to find a primitive cell constrained by these planes.
One issue here is that we rotated a big cube of atoms at a weird angle, and this means that, looking along our $[\bar{2}11], [111], [01\bar{1}]$ direction, some atoms will be hidden behind others in a non-periodic manner, leading to a deceptive arrangement of atoms. You'll have to delete some atoms to make a proper decision as to where the planes go.
I picked the following positions of the planes.
| $h$ | $k$ | $l$ | distance ($\times d$) |
| --- | --- | --- |:---------------------:|
| -2 | 1 | 1 | -4 |
| -2 | 1 | 1 | 2 |
| 1 | 1 | 1 | -1 |
| 1 | 1 | 1 | 2 |
| 0 | 1 | -1 | -1 |
| 0 | 1 | -1 | 1 |
> [!WARNING]
> There are clearly multiple possible primitive unit cells to build here (related by an origin shift). There are several possible atomic terminations to cleave the sides of the wire ($(111)$ and $(01\bar{1})$ surface), I just picked one of them.
> [!TIP]
> Depending on your application, it might be necessary to expose certain atoms while constructing the wire. However, since we're constructing a "nanowire primitive cell", you can easily shift the origin of that cell before building a large $N\times N\times 1$ supercell and adding vacuum.
>
> The shape of the cell in the $z$-direction does not really matter since it will be periodic in the final wire, but I tried to make it as nice as possible.
#### New Lattice Vectors
This is the primitive cell we were going for. We need to pick the lattice vectors of our new unit cell, and for that we need four corner atoms, an example is shown below.
![[Sr3PbO_primnw_trimmed 8.png]]
![[Sr3PbO_primnw_trimmed 7.png]]
The atoms are labelled, and their coordinates are shown below:
$
\begin{align}
\tilde{o} &= (1, -3/2, -1/2), \\
\tilde{a} &= (2, -1/2, 1/2),\\
\tilde{b} &= (1, -1/2, -3/2), \\
\tilde{c} &= (-1, -1/2, 1/2), \\
\end{align}
$
Thus, we get the following lattice vectors from the origin $\tilde{o}$
$
\begin{align}
\tilde{\mathbf{a}} &= (1, 1, 1), \\
\tilde{\mathbf{b}} &= (0, 1, -1), \\
\tilde{\mathbf{c}} &= (-2, 1, 1)
\end{align}
$
Since the unit cell is cubic, the transformation matrix is then simply:
$
T = \begin{pmatrix}
1 & 1 & 1 \\
0 & -1 & 1 \\
2 & -1 & -1
\end{pmatrix}
$
This is more or less equivalent to what we had before, you just need to multiply the bottom two rows by $-1$, which doesn't actually change the transformed cell. We could have got the exact same matrix as before by picking different atoms.
## Constructing the "primitive" nanowire cell
All we need to do now is transform the *primitive* unit cell using $T$ above.
> [!TIP]
> You need to take the transpose of the transformation matrix before using it in VESTA, it uses a different convention.
> [!WARNING]
> VESTA tends to have issues with transformation matrices with half integer entries (it generates too many duplicate atoms). `pymatgen` and `ase` can't handle non-integer transformation matrices.
> If you end up with such a matrix, multiply it by 2 and use that as your new transformation matrix, then trim the resultant cell by hand.
This results in the following structure, which is fully periodic and identical to the one we had trimmed.
![[Sr3PbO_primnw_pymatgen.png]]
>[!TIP]
> An easy way to check your nanowire supercell is passing it through [FINDSYM](https://stokes.byu.edu/iso/findsym.php). It should reproduce the original primitive unit cell.
# Constructing the Nanowire
## Manually using VESTA
If you don't want to automate the whole process, you can construct the nanowire using VESTA. Save the "primitive" nanowire cell in a new `.cif` or `POSCAR` file and open it in VESTA.
First, we create a supercell of the "primitive" nanowire cell, e.g., $5\times5\times1$ (this depends on the cross section you need). You can fiddle around with the terminations in the $a$ and $b$ directions if needed. Save this cell into a VASP `POSCAR` with *cartesian* coordinates.
The next step is a little counterintuitive. We need the wire to be periodic in the $c$ direction, so we have to delete the layer of atoms at $c = 1$, since they are just a periodic image of the layer at $c = 0$. Export the structure to `.xyz` and do not save the hidden atoms.
The `.xyz` will have all the atoms at the $a = 1$ and $b = 1$ planes, which are missing from the `POSCAR`. Open the cartesian `POSCAR` in a text editor and replace the list of atoms by the one in the `.xyz` file (make sure you delete the first column with atomic symbols after you paste the new list, or move it to be the fourth column after the $z$ coordinate).
Now, we can add some vacuum to the $a$ and $b$ direction, e.g., $15\, \text{Å}$ . Add $15$ to the $a$ and $b$ lattice constants in the `POSCAR` and we're almost done.
The final step is changing the number of atoms in the `POSCAR`. You can use a terminal command on the `.xyz` file to count the atoms in it, e.g.,
```bash
$ grep -o -i Sr my_struct.xyz | wc -l
> 497
$ grep -o -i Pb my_struct.xyz | wc -l
> 177
$ grep -o -i O my_struct.xyz | wc -l
> 166
```
Subtract 1 from these numbers (the first line in a `.xyz` file contains all the atomic symbols), and change the `POSCAR` to these numbers.
Finally, we get the following:
![[my_struct 1.png]]
The atoms on the $a = 1$ and $b = 1$ planes are just periodic images and do not actually show up in the calculation.
> [!SUMMARY]
> 1. Construct a supercell from the "primitive" nanowire cell using VESTA
> 2. Export to *cartesian* `POSCAR`.
> 3. Delete the top layer of atoms at $c = 1$.
> 4. Export to `.xyz` and do not save hidden atoms.
> 5. Copy the list of atoms form the `.xyz` file to the `POSCAR`, and get rid of the first column of atomic symbols.
> 6. Add the size of the vacuum to the lattice constants in the $a$ and $b$ directions.
> 7. Count the atoms of each species in the `.xyz` and adjust the `POSCAR`
## Automatically using `pymatgen`
My applications requires testing quite a few nanowire structures, and armed with the transformation matrix above, we can automate the whole process using `pymatgen`.
I will improve this script and post it on GitHub soon. I will include the link here later.
```python
import numpy as np
from pymatgen.io.vasp import Poscar
from pymatgen.io.cif import CifWriter
from pymatgen.core.structure import Structure
import math
from operator import itemgetter
from pymatgen.core.operations import SymmOp
import pandas as pd
# Size of vacuum in a and b directions
vac_a = 15
vac_b = 15
# Size of supercell in a and b directions
N_a = 5
N_b = 5
# Convert primitive cell to nanowire cell.
prim = Structure.from_file("Sr3PbO_primitive.cif")
T = [[1, 1, 1], [0, -1, 1], [2, -1, -1]]
prim.make_supercell(scaling_matrix=T)
# Shift Origin (Optional)
#origin_shift = [-1/6, 0, 1/3]
#symop = SymmOp.from_rotation_and_translation(rotation_matrix=np.eye(3), translation_vec=origin_shift)
#prim.apply_operation(symop, fractional=True)
#sg = SpacegroupAnalyzer(prim)
#prim = sg.get_refined_structure()
# For some reason, the rersultant NW cell doesn't have its axes aligned with cartesian x y z, so we fix that
# Figure out the transformation matrix
a = np.array([np.linalg.norm(T[0]), 0, 0])
b = np.array([0, np.linalg.norm(T[1]), 0])
c = np.array([0, 0, np.linalg.norm(T[2])])
final_vectors = np.stack((a, b, c))
new_T = np.matmul(np.linalg.inv(T), final_vectors)
# Apply Transformation
symop = SymmOp.from_rotation_and_translation(rotation_matrix=np.transpose(new_T))
prim.apply_operation(symop)
# Save "primitive" NW cell for checking
Poscar(prim).write_file('Sr3PbO_primnw_pymatgen.vasp', direct=True, significant_figures=10)
# Make the supercell to form the actual wire
prim.make_supercell([N_a,N_b,1])
# You can also save the N_a*N_b*1 supercell if you want to check stuff
#Poscar(prim).write_file('Sr3PbO_primnw_5x5_pymatgen.vasp', direct=True, significant_figures=10)
# Create a cartesian POSCAR string from the supercell, they're easy to work with
new_poscar = Poscar(prim)
poscar_string = new_poscar.get_string(direct=False, significant_figures=10)
string_array = poscar_string.split('\n')
n_atoms = len(string_array) - 8 # Assumes VASP 5 mode
string_array.pop() # Get rid of empty line
# Change lattice constants to add vacuum
a = np.array(string_array[2].split()).astype(float)
new_a = a + [vac_a, 0, 0]
b = np.array(string_array[3].split()).astype(float)
new_b = b + [0.0, vac_b, 0]
# Write them out to the array with the POSCAR lines.
string_array[2] = ' ' + np.array2string(new_a, formatter={'float_kind':lambda x: "%.10f" % x}, separator = " ")[1:-2]
string_array[3] = ' ' + np.array2string(new_b, formatter={'float_kind':lambda x: "%.10f" % x}, separator = " ")[1:-2]
# Extract atomic positions and shift them
# Atomic position in POSCAR start on the 9th line for VASP 5+
pos_string_array = string_array[8:]
# Symbols of the elements in the structure
el_sym = string_array[5].split()
# Counter for how many atoms we add of each type
extra_atoms = {key: 0 for key in el_sym}
# Shift and duplicate atoms
for n, pos_string in enumerate(pos_string_array):
pos = np.array(pos_string.split()[0:3]).astype(float) # atom position
atom_sym = ' ' + pos_string.split()[-1] # atom symbol
bool_0_0 = math.isclose(pos[0], 0, abs_tol=1e-3) # Does it lie on the a = 0 plane?
bool_0_1 = math.isclose(pos[0], a[0], abs_tol=1e-3) # Does it lie on the a = 1 plane?
bool_1_0 = math.isclose(pos[1], 0, abs_tol=1e-3) # Does it lie on the b = 0 plane?
bool_1_1 = math.isclose(pos[1], b[1], abs_tol=1e-3) # Does it lie on the b = 1 plane?
if bool_0_1: # it's on the a = 1 plane, need to shift it to a = 0
print(f"WARNING, ATOM AT WRONG SIDE, pos = {pos}, shifted to {pos - a}")
pos -= a
pos_str = ' ' + np.array2string(pos, formatter={'float_kind':lambda x: "%.10f" % x}, separator = " ")[1:-2] + atom_sym
string_array[n+8] = pos_str
bool_0_0 = True # Now it's on the a = 0 plane
if bool_1_1: # it's on the b = 1 plane, need to shift it to b = 0
print(f"WARNING, ATOM AT WRONG SIDE, pos = {pos}, shifted to {pos - b}")
pos -= b
pos_str = ' ' + np.array2string(pos, formatter={'float_kind':lambda x: "%.10f" % x}, separator = " ")[1:-2] + atom_sym
string_array[n+8] = pos_str
bool_1_0 = True # Now it's on the b = 0 plane
if bool_0_0 and (not bool_1_0): # This atom lies on the a = 0 plane but not on the hinge (a=0,b=0)
new_pos = pos + [a[0], 0, 0]
pos_str = ' ' + np.array2string(new_pos, formatter={'float_kind':lambda x: "%.10f" % x}, separator = " ")[1:-2] + atom_sym
string_array.append(pos_str)
extra_atoms[atom_sym.strip()] += 1 # Count the extra atom we added
elif bool_1_0 and (not bool_0_0): # This atom lies on the b = 0 plane but not on the hinge (a=0,b=0)
new_pos = pos + [0, b[1], 0]
pos_str = ' ' + np.array2string(new_pos, formatter={'float_kind':lambda x: "%.10f" % x}, separator = " ")[1:-2] + atom_sym
string_array.append(pos_str)
extra_atoms[atom_sym.strip()] += 1 # Count the extra atom we added
elif (bool_0_0) and (bool_1_0): # This atom is on the hinge, need 3 copies
# To a = 1, b = 0 hinge
new_pos = pos + [a[0], 0, 0]
pos_str = ' ' + np.array2string(new_pos, formatter={'float_kind':lambda x: "%.10f" % x}, separator = " ")[1:-2] + atom_sym
string_array.append(pos_str)
# To b = 1, a = 0 hinge
new_pos = pos + [0, b[1], 0]
pos_str = ' ' + np.array2string(new_pos, formatter={'float_kind':lambda x: "%.10f" % x}, separator = " ")[1:-2] + atom_sym
string_array.append(pos_str)
# To b = 1, a = 1 hinge
new_pos = pos + [a[0], b[1], 0]
pos_str = ' ' + np.array2string(new_pos, formatter={'float_kind':lambda x: "%.10f" % x}, separator = " ")[1:-2] + atom_sym
string_array.append(pos_str)
extra_atoms[atom_sym.strip()] += 3 # Count the extra 3 atoms we added
# OPTIONAL: Symmetrize the vacuum.
for n, pos_string in enumerate(string_array[8:]):
pos = np.array(pos_string.split()[0:3]).astype(float)
new_pos = pos + [vac_a/2, vac_b/2, 0]
atom_sym = ' ' + pos_string.split()[-1]
pos_str = ' ' + np.array2string(new_pos, formatter={'float_kind':lambda x: "%.10f" % x}, separator = " ")[1:-2] + atom_sym
string_array[n+8] = pos_str
# Adjust the number of atoms in the POSCAR to account for the periodic copies
num_atoms = np.array(string_array[6].split()).astype(float) + np.array(list(extra_atoms.values()))
string_array[6] = ' '.join(map(str, num_atoms.astype(int)))
# Now we need to sort the POSCAR so the atoms are in the expected order.
atoms = [mystring.split() for mystring in string_array[8:]]
df = pd.DataFrame(atoms)
df.columns = ["x", "y", "z", "sym"]
sorting_order = dict(zip(el_sym, list(range(0,len(el_sym)))))
df['order'] = df.sym.map(sorting_order)
sorted_df = df.sort_values('order').drop('order',axis=1)
new_string_array = string_array[0:8] + [' '.join(mylist) for mylist in sorted_df.to_numpy()]
# Convert to pymatgen structure
final_wire = Structure.from_str('\n'.join(new_string_array), fmt='poscar')
# Print in different formats
Poscar(final_wire).write_file('Sr3PbO_wire_final.vasp', direct=True, significant_figures=10)
CifWriter(final_wire).write_file('Sr3PbO_wire_final.cif')
```
The resulting structure is shown below:
![[Sr3PbO_wire_final.png]]
![[Sr3PbO_wire_final 1.png]]
![[Sr3PbO_wire_final 3.png]]