voronoi图代码_周期边界Voronoi在分子动力学MD中的应用（生成多晶石墨烯）

感谢最近一些在做MD和随机网络的朋友的关注和私信，我决定在这里陆续分享一些相关MD的小知识，原因是这些内容在国内的很难搜索到，资料非常匮乏，同时也希望我可以记录自己的学习历程。Voronoi的生成不是难点，程序和介绍都非常多，这篇文章的特点主要区别是运用了周期性边界条件（PBC），以便进行后期的MD仿真，最后附有演化视频和几个Voronoi源代码。

MD的编程应用对我们的发展至关重要，工业界CAD和CAE已经被老外牢牢掌握，经过了数10年的发展，我们很难在上面进行超越，如果想用这些软件去发表论文和工业设计，都要老老实实的给老外交钱，一般这些软件都价格不菲，最重要的是一旦国外停止供应，就一个matlab都足以给我们造成很大的打击。然而在物理化学领域，第一性原理和MD几乎各占半边天，也有现在逐渐相结合的aiMD，幸运的是这些软件都需要很大的改进，目前还都不完善，没有一个非常强大且垄断的商业软件出现，所以我们国家在这方面有着巨大的发展潜力。

好了，废话不多说了，本文说说随机的二维Voronoi图（附带周期边界条件）如何变成石墨烯，通过这样产生的石墨烯脱离了完美晶体模型的记忆性，更加利于进一步的MD仿真！Voronoi又称“泰森多边形”，在模拟计算领域可谓是大名鼎鼎，广泛用在气象，地理，交通，晶体生长领域，长相如下：

一个平面被分成了多个区域，而且每个区域中都仅包含一个随机点，这里要注意的是，可能会误解是先有的彩色区域而后有的黑点，实际上是现有的黑色随机点，然后进行德劳内三角划分得到了不同的彩色区域。这样的图形有一个性质就是每个区域内任何一点到该随机点的距离总是小于到其他区域内随机点的距离。

具体的介绍可以参见里面的高赞回答： 沃罗诺伊图（Voronoi Diagram，也称作Dirichlet tessellation，狄利克雷镶嵌 ）是怎样的？ 最好有图说明。

方法不再赘述，我这里的Voronoi使用了首先生成一堆随机点的坐标（x, y)，x是在0~Lx间生成，y是在0~Ly之间生成，所以一定要记下这里的Lx和Ly，这是运用到周期边界中的关键因素，而且以便后期盒子松弛使用。

下面产生100个随机点，为了运用PBC，实际需要400个点，因为要把每个方位的点都按周期Lx和Ly移动到相应的位置上，例如：

上图中红色虚线中的是我们的仿真盒子，它在x和y方向的长度分别为Lx和Ly，在虚线的区域中放入N个原子，为了使用PBC，我们需要再加3N个原子按如图所示放在其他各个方向，从而构成了具有PBC的散点图，如下：

红色实线为我们的仿真box，里面100个原子，外面的为周期性复制，接下来进行Delaunay三角划分，将散点图构成三角形，然后找出其外接圆的圆心，此处不再赘述，然后将外接圆的圆心都链接起来，就得到了泰森多边形，每个区域内仅仅包含一个散点。

每个外接圆圆心都只与三个点相连，正好满足石墨烯中每个C原子外面有3个分子键。然后将随机点去掉，只留下box中的Voronoi图：

之前box内的100个散点图这里变成了200个C原子，然后给每个原子标号，记录下每个键两边连接的原子信息，还有每个原子相邻的3个原子信息，放入我们的MD程序中进行能量弛豫，这样就可以得到原始的石墨烯样本：

如果扩展到三维，只需要随机改变部分c原子的z坐标：

弯弯曲曲的看起来，不过这只是最初的随机模型，如果进一步生成多晶石墨烯还需要随机换键处理。

下面是二维voronoi弛豫的视频：

知乎视频​www.zhihu.com

附上一个简单的Voronoi图的python代码：

from PIL import Image
import random
import math

def generate_voronoi_diagram(width, height, num_cells):
	image = Image.new("RGB", (width, height))
	putpixel = image.putpixel
	imgx, imgy = image.size
	nx = []
	ny = []
	nr = []
	ng = []
	nb = []
	for i in range(num_cells):
		nx.append(random.randrange(imgx))
		ny.append(random.randrange(imgy))
		nr.append(random.randrange(256))
		ng.append(random.randrange(256))
		nb.append(random.randrange(256))
	for y in range(imgy):
		for x in range(imgx):
			dmin = math.hypot(imgx-1, imgy-1)
			j = -1
			for i in range(num_cells):
				d = math.hypot(nx[i]-x, ny[i]-y)
				if d < dmin:
					dmin = d
					j = i
			putpixel((x, y), (nr[j], ng[j], nb[j]))
	image.save("VoronoiDiagram.png", "PNG")
	image.show()
 
generate_voronoi_diagram(500, 500, 25)