Bug of ase.build.surface?

ASE
1

Hi,

When building fcc(111) surfaces, the surface module sometimes generates wrong c-axis length. The following test code aims to generate 3 (111) layers. However, for some input fcc lattice, only 2 layers are generated (with the 3rd coinciding with the 1st layer), and hence interlayer spacing and c-axis length are wrong. The code is “““

import numpy as np
from ase import Atoms
from ase.build import bulk
from ase.build import surface

def analyze_slab(a, n_layers=3):
“”“Analyze slab created from FCC with lattice constant a”“”

# Create bulk FCC
if 1:
    cu3au = Atoms('Cu3Au', 
                  scaled_positions=[(0.5, 0.5, 0), (0.5, 0, 0.5), 
                                   (0, 0.5, 0.5), (0, 0, 0)],
                  cell=[a, a, a], 
                  pbc=True)
#cu3au = bulk('Cu', 'fcc', a=a)  # test primitive 1-atom FCC
# Create (111) surface
atoms = surface(cu3au, (1, 1, 1), n_layers, vacuum=0, periodic=True)

# Calculate important parameters
cell = atoms.cell
interlayer_spacing = cell[2, 2] / n_layers
in_plane_area = abs(np.cross(cell[0, :2], cell[1, :2]))

print(f"Input fcc a = {a:.6f} ")
print(f"  c-axis length: {cell[2, 2]:.6f} ")
print(f"  Number of layers: {n_layers}")
print(f"  Interlayer spacing: {interlayer_spacing:.6f} ")
print(f"  Expected interlayer spacing: {a/np.sqrt(3):.6f} ")


return atoms

for a in [3.73, 3.74,3.75,3.76,3.77,3.78,3.79,3.8]:
#print(f"Case: a = {a:.4f} Å")
atoms1 = analyze_slab(a, n_layers=3)
print()

“““. The result is:

Input fcc a = 3.730000
c-axis length: 4.307033
Number of layers: 3
Interlayer spacing: 1.435678
Expected interlayer spacing: 2.153517

Input fcc a = 3.740000
c-axis length: 6.477870
Number of layers: 3
Interlayer spacing: 2.159290
Expected interlayer spacing: 2.159290

Input fcc a = 3.750000
c-axis length: 4.330127
Number of layers: 3
Interlayer spacing: 1.443376
Expected interlayer spacing: 2.165064

Input fcc a = 3.760000
c-axis length: 6.512511
Number of layers: 3
Interlayer spacing: 2.170837
Expected interlayer spacing: 2.170837

Input fcc a = 3.770000
c-axis length: 4.353221
Number of layers: 3
Interlayer spacing: 1.451074
Expected interlayer spacing: 2.176611

Input fcc a = 3.780000
c-axis length: 4.364768
Number of layers: 3
Interlayer spacing: 1.454923
Expected interlayer spacing: 2.182384

Input fcc a = 3.790000
c-axis length: 4.376315
Number of layers: 3
Interlayer spacing: 1.458772
Expected interlayer spacing: 2.188158

Input fcc a = 3.800000
c-axis length: 6.581793
Number of layers: 3
Interlayer spacing: 2.193931
Expected interlayer spacing: 2.193931

2

Hi Scidog,

I am not an ASE developer, but an avid user, so I’ll try my best to give a satisfying answer.

tl;dr: remove vacuum=0 from the argument list for a quick fix.

I think this is caused by a combination of setting vacuum=0 when you call surface and the subsequent atoms.wrap().

The surface is actually constructed correctly by surface in the beginning. However, setting vacuum=0 causes ASE to resize the lattice of your surface slab after it has been constructed to the exact bounding box of the atoms in the cell (by calling Atoms.center() internally). At this point, your cell is already wrong, since the bottom and top layer will overlap when you replicate the cell. Finally, when you call wrap, floating-point precision will cause some of the atoms (which sit right at the top of the cell) to be mapped to the bottom layer of the cell instead.

I snapped a couple of pictures to demonstrate what I mean:

image
image740×297 38.5 KB

The easy fix for your setup is to remove the vacuum=0 from the argument list:

image
image374×326 1.71 KB

However, when I looked at this I realized that the same problems cause an even bigger issue down the line when vacuum>0. Take a look at this slab with vacuum=5, which is not what we want:

image
image272×408 1.3 KB

As a fix for this, one can truncate the precision errors by changing ASE’s build tool like this:

--- a/ase/build/general_surface.py
+++ b/ase/build/general_surface.py
@@ -76,6 +76,9 @@ def build(lattice, basis, layers, tol, periodic):
     surf = lattice.copy()
     scaled = solve(basis.T, surf.get_scaled_positions().T).T
     scaled -= np.floor(scaled + tol)
+
+    scaled[abs(scaled) < tol] = 0.0
+
     surf.set_scaled_positions(scaled)
     surf.set_cell(np.dot(basis, surf.cell), scale_atoms=True)
     surf *= (1, 1, layers)

In principle, I could open up a merge request with this change. But since I am not an ASE developer, I would appreciate a second opinion :blush:

Best
Alex