1.LAMMPS介绍

• LAMMPS(Large-scale Atomic/Molecular Massively Parallel Simulator)，主要用于分子动力学模拟，但也可用于能量最小化和蒙特卡洛模拟，可模拟几个粒子到几十亿粒子的体系。
• 代码开源免费，使用C++编写，目前由美国桑迪亚国家实验室(Sandia National Laboratory，SNL)进行维护。
• 可在笔记本/台式机上运行，可安装在MacOS、Linux和Windows上。
• 常见的输入文件：
  1.structure.dat(data文件，包含体系大小、原子数、原子类型、各类型质量、各原子位置，也可包含原子速度)
  2.filename.in(input文件，包含运行程序的参数，也可包含势函数或经验势)
  3.pot.mod(potential文件，包含有关势参数)
• 常见的输出文件：
  1.final-structure.dat(data文件：与输入data文件相同，但包含最终体系的信息。可以作为另一个模拟的初始输入data文件)
  2.md.lammpstrj(轨迹文件：每N步的原子位置，用于进一步可视化分析)
  3.log.lammps(日志文件：存储所有LAMMPS执行的操作。当出现错误时，可以查看)
  4.out.filename(输出文件：包含所执行的操作的摘要，错误信息，以及热力学数据)

2.使用python编写data文件

1.使用python编写简单data文件

#使用python创建data文件
import numpy as npn_atoms = 1000                          #创建原子数
system_size = 20.0                      #设置体系长度，埃米，假设cubic体系positions = []                          #生成原子坐标
for i in range(n_atoms):                #随机生成，非晶态positions.append(np.random.rand(3)*system_size)with open("random.data","w") as fdata:                 #编写lammps的data文件fdata.write("random atoms - written\n\n")                #首行为注释行#头部fdata.write("{} atoms\n".format(n_atoms))                #指定原子数fdata.write("{} atom types\n".format(1))                   #指定原子类型数fdata.write("{} {} xlo xhi\n".format(0.0,system_size))              #指定盒子尺寸fdata.write("{} {} ylo yhi\n".format(0.0,system_size))fdata.write("{} {} zlo zhi\n".format(0.0,system_size))fdata.write("\n")#原子部分fdata.write("Atoms\n\n")for i,pos in enumerate(positions):                        #写入坐标fdata.write("{} 1 {} {} {}\n".format(i+1,*pos))

2.使用python编写晶体data文件

#使用python创建data文件
import numpy as nplattice_parameter = 4.046                 #晶格常数，埃米
basis = np.array([[1.0,0.0,0.0],            #基础FCC立方[0.0,1.0,0.0],[0.0,0.0,1.0]])*lattice_parameter
base_atoms = np.array([[0.0,0.0,0.0],[0.5,0.5,0.0],[0.5,0.0,0.5],[0.0,0.5,0.5]])*lattice_parameter
system_size = 10                    #设置体系长度，晶格单位，假设cubic体系positions = []                      #生成原子坐标
for i in range(system_size):for j in range(system_size):for k in range(system_size):base_position = np.array([i,j,k])cart_position = np.inner(basis.T,base_position)        # np.inner()计算内积for atom in base_atoms:positions.append(cart_position + atom)with open("crystalline.data","w") as fdata:            #编写lammps的data文件fdata.write("crystalline Al atoms - written\n\n")              #首行为注释行#头部fdata.write("{} atoms\n".format(len(positions)))            #指定原子数fdata.write("{} atom types\n".format(1))              #指定原子类型数fdata.write("{} {} xlo xhi\n".format(0.0,system_size*lattice_parameter))              #指定盒子尺寸fdata.write("{} {} ylo yhi\n".format(0.0,system_size*lattice_parameter))fdata.write("{} {} zlo zhi\n".format(0.0,system_size*lattice_parameter))fdata.write("\n")#原子部分fdata.write("Atoms\n\n")for i,pos in enumerate(positions):                    #写入坐标fdata.write("{} 1 {} {} {}\n".format(i+1,*pos))

3.使用python编写晶体(旋转)data文件

#使用python创建data文件
import numpy as np
import itertools as its
import math#新代码，根据轴和角度旋转矩阵
def rotation_matrix(axis,angle):a = math.cos(angle/2.0)s = math.sin(angle/2.0)b = axis[0]*sc = axis[1]*sd = axis[2]*sreturn np.array([[a*a+b*b-c*c-d*d,2*(b*c-a*d),2*(b*d+a*c)],[2*(b*c+a*d),a*a+c*c-b*b-d*d,2*(c*d-a*b)],[2*(b*d-a*c),2*(c*d+a*b),a*a+d*d-b*b-c*c]])lattice_parameter = 4.046                  #晶格常数，埃米
basis = np.array([[1.0,0.0,0.0],           #基础FCC立方[0.0,1.0,0.0],[0.0,0.0,1.0]])*lattice_parameter
base_atoms = np.array([[0.0,0.0,0.0],[0.5,0.5,0.0],[0.5,0.0,0.5],[0.0,0.5,0.5]])*lattice_parameter
#旋转后的基础FCC立方
axis = [1,1,0]             #轴
angle = np.pi/6            #角度
rot = rotation_matrix(axis,angle)
basis = np.matmul(rot.T,basis)                    # np.matmul()矩阵相乘
base_atoms = np.array([np.matmul(rot,atom) for atom in base_atoms])system_size = 20                  #设置体系长度，晶格埃米，假设cubic体系
system_basis = np.array([[1.0,0.0,0.0],[0.0,1.0,0.0],[0.0,0.0,1.0]])*system_size
inv_basis = np.linalg.inv(basis)                          # np.linalg.inv()计算逆协方差矩阵
inv_system = np.linalg.inv(system_basis)
corners = np.array(list(its.product([0.0,system_size],repeat = 3)))          # its.product() 迭代循环 3次
extents = np.array([np.matmul(inv_basis.T,corner) for corner in corners])
mins = [int(np.min([m,0])) for m in np.floor(np.min(extents,axis = 0))]             # np.floor() 向下取整
maxs = [int(np.max([m,0]))+1 for m in np.ceil(np.max(extents,axis = 0))]             # np.ceil() 向上取整positions = []                    #生成原子坐标
for i in range(mins[0],maxs[0]):for j in range(mins[1],maxs[1]):for k in range(mins[2],maxs[2]):base_position = np.array([i,j,k])cart_position = np.inner(basis.T,base_position)          # np.inner()计算内积for atom in base_atoms:atom_pos = cart_position + atomsystem_pos = np.matmul(inv_system.T,atom_pos)if all(system_pos >= 0.0) and all(system_pos < 1.0):        #如果该原子在盒子中positions.append(atom_pos)with open("rotate_crystalline.data","w") as fdata:        #编写lammps的data文件fdata.write("rotate crystalline Al atoms - written\n\n")              #首行为注释行#头部fdata.write("{} atoms\n".format(len(positions)))            #指定原子数fdata.write("{} atom types\n".format(1))                   #指定原子类型数fdata.write("{} {} xlo xhi\n".format(0.0,system_size))              #指定盒子尺寸fdata.write("{} {} ylo yhi\n".format(0.0,system_size))fdata.write("{} {} zlo zhi\n".format(0.0,system_size))fdata.write("\n")#原子部分fdata.write("Atoms\n\n")for i,pos in enumerate(positions):                    #写入坐标fdata.write("{} 1 {} {} {}\n".format(i+1,*pos))