如何找到扭曲双层石墨烯超晶格的晶胞？

Question

我试图编写一个代码，在其中生成扭曲双层石墨烯的坐标，然后找到两层之间的匹配点，以获得最终的莫尔图案。但是有些东西在我的代码中不起作用，它没有给我我需要的解决方案。我在这里找不到错误。

import numpy as np
import matplotlib.pyplot as plt
from itertools import groupby 
import math

d0 = 0.3330630630630631 # interlayer distance
a0 = 0.15469469469469468 # distance between the two basis atoms of graphene

theta = 15*np.pi/180 # the twist angle between both layers in degree

nmax = 10 

# lattice vectors sublattice 1
a1= np.array([3/2*a0,3**0.5/2*a0,0])
a2= np.array([3/2*a0,-3**0.5/2*a0,0])

# lattice vectors sublattice 2
b1 = np.array([ np.cos(theta) * a1[0] - np.sin(theta) * a1[1], np.sin(theta) * a1[0] + np.cos(theta) * a1[1], 0. ])
b2 = np.array([ np.cos(theta) * a2[0] - np.sin(theta) * a2[1], np.sin(theta) * a2[0] + np.cos(theta) * a2[1], 0. ])


#Coordinates sublattice 1 & 2 for layers 1 & 2
coords1a = np.array([i * a1 + j * a2 for i in range(-nmax,nmax+1) for j in range(-nmax,nmax+1)])
coords1b = np.array([i * a1 + j * a2 + [a0,0.,0.] for i in range(-nmax, nmax+1) for j in range(-nmax, nmax+1)])


coords2a = np.array([i * b1 + j * b2 + [0.,0.,d0] for i in range(-nmax,nmax+1) for j in range(-nmax,nmax+1)])
coords2b = np.array([i * b1 + j * b2 + [np.cos(theta)*a0, np.sin(theta)*a0,d0] for i in range(-nmax,nmax+1) for j in range(-nmax,nmax+1)])

coords1 = np.concatenate((coords1a, coords1b), axis=0)
coords2 = np.concatenate((coords2a, coords2b), axis=0)


# to find the matching coordinates of the two layers that are supposed to form the moire pattern

for a in range(len(coords1)):
    for b in range(len(coords2)):

        dis1 = coords1[a][0] - coords2[b][0]
        dis2 = coords1[a][1] - coords2[b][1]
        distance1 = np.sqrt(coords1[a][0]**2 + coords1[a][1]**2)
        distance2 = np.sqrt(coords2[b][0]**2 + coords2[b][1]**2)
        if abs(dis1) < small_variable:
            if abs(dis2) < small_variable:
                if distance1 == distance2: 
                    matching_coords.append(coords1[a])
                    print(coords1[a], ' is matching with ', coords2[b])

我尝试更改向量并使用小变量的不同值，但到目前为止没有任何效果。