Tutorial - MPB DocumentationMPB官网给出了一个示例程序说明最大化三角晶格光子晶体第一带隙的方法代码是import mathimport meep as mpfrom meep import mpbimport numpy as npimport matplotlib.pyplot as pltnum_bands 2resolution 32geometry_lattice mp.Lattice(sizemp.Vector3(1, 1),basis1mp.Vector3(math.sqrt(3)/2, 0.5),basis2mp.Vector3(math.sqrt(3)/2, -0.5))geometry [mp.Cylinder(0.2, materialmp.Medium(epsilon12))]k_points [mp.Vector3(), # Gammamp.Vector3(y0.5), # Mmp.Vector3(-1./3, 1./3), # Kmp.Vector3(), # Gamma]k_points mp.interpolate(4, k_points)ms mpb.ModeSolver(geometrygeometry,geometry_latticegeometry_lattice,k_pointsk_points,resolutionresolution,num_bandsnum_bands)ms.mesh_size 7def first_tm_gap(r):ms.geometry [mp.Cylinder(r, materialmp.Medium(epsilon12))]ms.run_tm()return -1 * ms.retrieve_gap(1) # return the gap from TM band 1 to TM band 2from scipy.optimize import minimize_scalarresult minimize_scalar(first_tm_gap, methodbounded, bounds[0.1, 0.5], options{xatol: 0.1})print(radius at maximum: {}.format(result.x))print(gap size at maximum: {}.format(result.fun * -1))按理来说运行完应该出现radius at maximum: 0.176393202250021 gap size at maximum: 48.6252611051049但是我运行以后显示radius at maximum: 0.15835921350012622gap size at maximum: 48.15177339844839和官网的结果不同我的代码完全是从网站上复制下来的我思考了一下应该是“options{xatol: 0.1}”设置的不对容差率设置的太大了改成options{xatol: 0.01})就对了能出现正确结果了这说明官网的程序未必都正确我们接下来看呈现点缺陷超晶胞态的例程按照官网代码原原本本复制以后import mathimport meep as mpfrom meep import mpbimport numpy as npimport matplotlib.pyplot as pltnum_bands 50resolution 16ms.geometry_lattice mp.Lattice(sizemp.Vector3(5, 5))ms.geometry [mp.Cylinder(0.2, materialmp.Medium(epsilon12))]ms.geometry mp.geometric_objects_lattice_duplicates(ms.geometry_lattice, ms.geometry)ms.geometry.append(mp.Cylinder(0.2, materialmp.air))ms.k_points [mp.Vector3(0.5, 0.5)]ms mpb.ModeSolver(geometryms.geometry,geometry_latticems.geometry_lattice,k_pointsms.k_points,resolutionresolution,num_bandsnum_bands)ms.run_tm()ms.get_dfield(25) # compute the D field for band 25ms.compute_field_energy() # compute the energy density from Dc mp.Cylinder(1.0, materialmp.air)print(energy in cylinder: {}.format(ms.compute_energy_in_objects([c])))报错了原来是所有的ms.都应该被去除。改成import mathimport meep as mpfrom meep import mpbimport numpy as npimport matplotlib.pyplot as pltnum_bands 50resolution 16geometry_lattice mp.Lattice(sizemp.Vector3(5, 5))geometry [mp.Cylinder(0.2, materialmp.Medium(epsilon12))]geometry mp.geometric_objects_lattice_duplicates( geometry_lattice, geometry)geometry.append(mp.Cylinder(0.2, materialmp.air))k_points [mp.Vector3(0.5, 0.5)]ms mpb.ModeSolver(geometry geometry,geometry_lattice geometry_lattice,k_points k_points,resolutionresolution,num_bandsnum_bands)run_tm()get_dfield(25) # compute the D field for band 25compute_field_energy() # compute the energy density from Dc mp.Cylinder(1.0, materialmp.air)print(energy in cylinder: {}.format( compute_energy_in_objects([c])))以后就可以得到和官网一样的结果了。再看看 使用算法调节点缺陷的特征频率位置官网教程是import mathimport meep as mpfrom meep import mpbimport numpy as npimport matplotlib.pyplot as pltnum_bands 1resolution 16geometry_lattice mp.Lattice(sizemp.Vector3(5, 5))geometry [mp.Cylinder(0.2, materialmp.Medium(epsilon12))]geometry mp.geometric_objects_lattice_duplicates( geometry_lattice, geometry)geometry.append(mp.Cylinder(0.2, materialmp.air))k_points [mp.Vector3(0.5, 0.5)]ms mpb.ModeSolver(geometry geometry,geometry_lattice geometry_lattice,k_points k_points,resolutionresolution,num_bandsnum_bands)ms.target_freq (0.2812 0.4174) / 2ms.tolerance 1e-8from scipy.optimize import ridderold_geometry ms.geometry # save the 5x5 grid with a missing roddef rootfun(eps):# add the cylinder of epsilon eps to the old geometry:ms.geometry old_geometry [mp.Cylinder(0.2, materialmp.Medium(epsiloneps))]ms.run_tm() # solve for the mode (using the targeted solver)print(epsilon {} gives freq. {}.format(eps, ms.get_freqs()[0]))return ms.get_freqs()[0] - 0.314159 # return 1st band freq. - 0.314159rooteps ridder(rootfun, 1, 12)print(root (value of epsilon) is at: {}.format(rooteps))官网结果是我运行出来的倒数第二行是psilon 5.288831274341146 gives freq. 0.314158999981831最优的介电常数和官网不一致5.28和5.41差的太多了我现在还不知道是什么原因导致的最后看看查看特定模式电磁波分布的代码import mathimport meep as mpfrom meep import mpbimport matplotlib.pyplot as pltimport numpy as np #我发现这一行官网没有需要自己加上num_bands 8resolution 32geometry_lattice mp.Lattice(sizemp.Vector3(1, 1),basis1mp.Vector3(math.sqrt(3)/2, 0.5),basis2mp.Vector3(math.sqrt(3)/2, -0.5))geometry [mp.Cylinder(0.2, materialmp.Medium(epsilon12))]k_points [mp.Vector3(), # Gammamp.Vector3(y0.5), # Mmp.Vector3(-1./3, 1./3), # Kmp.Vector3(), # Gamma]k_points mp.interpolate(4, k_points)ms mpb.ModeSolver(geometrygeometry,geometry_latticegeometry_lattice,k_pointsk_points,resolutionresolution,num_bandsnum_bands)efields []def get_efields(ms, band):efields.append(ms.get_efield(band, bloch_phaseTrue))ms.run_tm(mpb.output_at_kpoint(mp.Vector3(1 / -3, 1 / 3), mpb.fix_efield_phase,get_efields))# Create an MPBData instance to transform the efieldsmd mpb.MPBData(rectifyTrue, resolution32, periods3)converted []for f in efields:# Get just the z component of the efieldsf f[..., 0, 2]converted.append(md.convert(f))for i, f in enumerate(converted):plt.subplot(331 i)plt.contour(converted_eps.T, cmapbinary) #官网上是带有.T的。但是我发现需要去掉.T来获得正确的图片不然的话得到的电磁波分布和EPSILON分布是旋转了90度的。plt.imshow(np.real(f).T, interpolationspline36, cmapRdBu, alpha0.9)plt.axis(off)plt.show()官网的结果是我得到的结果是虽然颜色表示有些不同官网的红色在我这是蓝色反之亦然但是内核是一样的。但是非常奇怪的是我的第78条能带的分布在官网结果中正号颠倒了不知道原因是啥。