I'm trying to find the adsorption binding sites of a lattice provided by a quantum espresso input file which has the following cell vectors:
kbasis=([[-5.14439960600149, 4.455181658787211, 0.0], [-3.858299704501117, -6.682772488180817, 0.0], [0.0, 0.0, -30.816597639823016]])
and the following atomic positions:
positions=([[ -3.8583 , -2.227591 , -12.00001486],
[ -3.86407716, -3.71598711, -14.09073627],
[ -5.1444 , -2.97679223, -16.18148849],
[ -3.8583 , -2.227591 , -18.2722099 ],
[ -6.4305 , 2.227591 , -12.00001486],
[ -6.43627716, 0.73919489, -14.09073627],
[ -2.5722 , -2.97679223, -16.18148849],
[ -6.4305 , 2.227591 , -18.2722099 ],
[ -5.1444 , 0. , -12.00001486],
[ -5.15017716, -1.48839611, -14.09073627],
[ -6.4305 , -0.74920123, -16.18148849],
[ -5.1444 , 0. , -18.2722099 ],
[ -2.5722 , 0. , -12.00001486],
[ -2.57797716, -1.48839611, -14.09073627],
[ -3.8583 , -0.74920123, -16.18148849],
[ -2.5722 , 0. , -18.2722099 ],
[ -9.0027 , -2.227591 , -12.00001486],
[ -5.15017716, 2.96678589, -14.09073627],
[ -1.2861 , -0.74920123, -16.18148849],
[ -9.0027 , -2.227591 , -18.2722099 ],
[ -3.8583 , 2.227591 , -12.00001486],
[ -3.86407716, 0.73919489, -14.09073627],
[ -5.1444 , 1.47838977, -16.18148849],
[ -3.8583 , 2.227591 , -18.2722099 ],
[ -5.1444 , -4.455182 , -12.00001486],
[ -1.29187716, 0.73919489, -14.09073627],
[ -2.5722 , 1.47838977, -16.18148849],
[ -5.1444 , -4.455182 , -18.2722099 ],
[ -7.7166 , 0. , -12.00001486],
[ -7.72237716, -1.48839611, -14.09073627],
[ -3.8583 , -5.20438323, -16.18148849],
[ -7.7166 , 0. , -18.2722099 ],
[ -6.4305 , -2.227591 , -12.00001486],
[ -6.43627716, -3.71598711, -14.09073627],
[ -7.7166 , -2.97679223, -16.18148849],
[ -6.4305 , -2.227591 , -18.2722099 ]])
The code I'm using is the following:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from pymatgen.analysis.adsorption import AdsorbateSiteFinder
from pymatgen.symmetry import analyzer
from pymatgen.symmetry.analyzer import SpacegroupOperations
from pymatgen.symmetry.analyzer import SpacegroupAnalyzer
from pymatgen import core
from pymatgen.core import Lattice
from pymatgen.core import Structure
from pymatgen.core import surface
import ase
from ase.io import read
tol=0.5
# Data extraction
data = read('enter image description herepwscf.in', format='espresso-in')
positions = data.get_positions()
labels_tmp = data.arrays['numbers']
labels=[];
for i in range(len(labels_tmp)):
if labels_tmp[i]==6:
labels=labels+["C"]
elif labels_tmp[i]==8:
labels=labels+["O"]
elif labels_tmp[i]==29:
labels=labels+["Cu"]
elif labels_tmp[i]==47:
labels=labels+["Ag"]
elif labels_tmp[i]==79:
labels=labels+["Au"]
kbasis = data.cell
# Lattice creation
surf_lattice=Lattice([[kbasis[0]],[kbasis[1]],[kbasis[2]]])
surf_structure=Structure(surf_lattice,labels,positions)
symm_finder=SpacegroupAnalyzer(surf_structure, symprec=0.1)
symm_ops = symm_finder.get_symmetry_operations(cartesian=True)
# Finding of the adsorption sites
ads_finder=AdsorbateSiteFinder(surf_structure)
ads=ads_finder.find_adsorption_sites(distance=1.8)
But the problem is that the binding sites I get are different from the correct ones, indeed if I compare, for example, the 'ontop' sites
ads_sites=np.array(ads['ontop'])
with the surface atomic positions
surf_positions=positions[positions[:,2]>max(positions[:,2])-tol]
I observe that these two set of points do not coincide!
plt.scatter(surf_positions[:,0],surf_positions[:,1],c='red',s=100)
plt.scatter(ads_sites[:,0],ads_sites[:,1],c='blue',s=50)
plt.show()
(see the image).
Please, help me, I have no idea on how to proceed.
Thank you
Axel
The positions of the possible adsorption sites are not correct.
