Bug of ase.build.surface?

[Scidog]

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

[aknoll]

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:

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

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:

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

Best
Alex

[askhl]

Hi,

Thanks for this extraordinary investigation @aknoll .

I had a brief look at the code, and we should probably apply this fix. However the code overall strikes me as a rather unsafe – if this fix is needed, I cannot really tell from only reading it that there aren’t other issues.

Another questionable thing is that the docstring claims that the builder method only works for “the conventional cell”, not the primitive one. I suppose that’s because the crystal directions 111 etc. are defined with respect to the conventional cell, but actually, since the code always does something specific, it should be possible to understand what that is and describe that. After all, it appears to do something useful also with non-conventional cells

But taking a step back: The code is based on layers, and in order for it to be reliable, it probably needs to know what layer each atom belongs to. If that is well defined, it can do the wrapping correctly. Otherwise it becomes only tolerance based, and there might be cases that slip through. Maybe the difficulty lies in the fact that multiple layers can come from each repetition of the basis, and some of those layers may need wrapping when producing the full system from the lattice repetitions. But it didn’t write which atoms went into the layer as they were generated.

[askhl]

By the way to be clear, a merge request would be highly welcome, particularly with a unittest which would fail before the fix and pass after the fix.

[aknoll]

Thank you @askhl, I agree with you that the code is a little bit obscure As a first step, I will open up a MR for this. If I can, I will add a unit test and extend the docstring, too.