Python-分子动力学模拟-Lammps-离子液体-统计受限体系粒子的个数

小弟我第一次写博客，课题组沿用十几年的fortran源代码，一个人一个声明，根本看不懂，网上跟着黑马学习了python，写了几个程序，天下的语言都好，我一直都觉得能解决问题的语言就是好语言，管它是白猫还是黑猫，python的优点在于读取你想要的数据这点比较方便，缺点就是效率低下，反正我是工作站，有的是时间等他处理数据。
言归正传，一般获取体系的粒子数是建模的第一步，不管你是NVT,NPT获取的也罢，都要统计受限区域的粒子个数，此博客的离子液体是[BMIM][BF4]和[BMIM][PF6],阳离子的中心原子采用的是N，阴离子的中心原子分别采用的是B和P。
可以交流学习一下，目前研二在读，一篇EA期刊IF=6.215，3月份接收的，
DOI: 10.1016/j.electacta.2020.135981
现在闲下来了，整理一下源代码。

图片

获取受限体系的平均粒子数

源码如下

注释的比较详细：
此源代码针对lammps输出的轨迹文件

# -*- coding: utf-8 -*-
# @Time : 2020/7/8 10:35 
# @Author :Jiabao_Zhu
# @File : read_confinedMolNumber.py 
# @Software: PyCharm
import pandas as pd
from numba import jit

LINES = 30000  # 看看是不是包含石墨烯的碳原子(不包含有30000) 33840
H_LINE = 9  # 表头,轨迹文件有几行空行
frames = 2002  # frame个数+1 因为range()函数左闭右开
rangeXfrom = -21.472  # XYZ的起始范围
rangeXto = 21.472
rangeYfrom = -19.905
rangeYto = 19.905
rangeZfrom = -18.375
rangeZto = 18.375

skiprows = [j * (LINES + H_LINE) + i for j in range(frames) for i in range(H_LINE)]
path = "./dumpNVT.lammpstrj"
# 取出数据,第2,3,4,5,6列,分别对应轨迹文件里的mol,type,x,y,z列
data = pd.read_csv(path, sep=' ', skiprows=skiprows, usecols=[1, 2, 3, 4, 5], header=None).values
f = open('confinedMolNumber.txt', 'w+', encoding='utf-8')  # 写入文件


@jit  # numba加速
def isInConfined_Anion(data):
    B_Number_List = []  # 创建符合条件的B原子个数的空列表,用来储存统计之后每一帧的个数
    for i in range(int(data.shape[0]/LINES)):  # data.shape[0]/LINES 指所有的数据行/一帧的行数,求出来就是多少帧的循环
        # 取出N1-bmim type 2 B-BF4 type 9 质心
        data_slice = data[i*LINES:(i+1)*LINES]  # 把所有帧的数据data进行切片,取出每一帧的数据进行处理
        atom_B = data_slice[data_slice[:, 1] == 9, 2:]  # data_slice[:,0]==9取出type为9的坐标
        B = atom_B
        B_index = []
        for j in range(B.shape[0]):  # shape[0]代表行数,shape[1]代表列数 B.shape[0]指取出来的B原子的行数=B的分子数
            # B-BF4边界判定
            if rangeXfrom < B[j][0] < rangeXto and rangeYfrom < B[j][1] < rangeYto and rangeZfrom < B[j][2] < rangeZto:
                B_index.append(j)
        B_Number = 0  # 符合条件的B原子个数的初始化
        for k in B_index:  # 遍历B_index列表,查看有多少个,
            B_Number += 1  # 一帧下的符合条件的个数
        # print(B_Number)
        B_Number_List.append(B_Number)  # 所有帧下的符合条件的个数添加到B_Number_List,如果有2帧，就应该有两个数据，2000帧，就应该有2000个数据
    # print(B_Number_List)
    B_Number_ = 0
    for l in B_Number_List:
        B_Number_ += l
    avg_B = str(B_Number_/(frames-1))  # 转化为str f.writelines(list_anion)只能识别str
    # print("阴离子平均个数为：", avg_B)
    list_anion = ["阴离子平均个数为：", avg_B]
    f.writelines(list_anion)


@jit
def isInConfined_Cation(data):
    N_Number_List = []
    for i in range(int(data.shape[0]/LINES)):
        # 取出N1-bmim type 2 B-BF4 type 9 质心
        data_slice = data[i*LINES:(i+1)*LINES]
        atom_N = data_slice[data_slice[:, 1] == 2, 2:]
        N = atom_N
        N_index = []
        for j in range(N.shape[0]):
            if rangeXfrom < N[j][0] < rangeXto and rangeYfrom < N[j][1] < rangeYto and rangeZfrom < N[j][2] < rangeZto:
                N_index.append(j)
        N_Number = 0
        for k in N_index:
            N_Number += 1
        N_Number_List.append(N_Number)
    N_Number_ = 0
    for l in N_Number_List:
        N_Number_ += l
    avg_N = str(N_Number_/((frames-1)*2))  # 转化为str
    # print("阳离子平均个数为：", avg_N)
    list_cation = ["阳离子平均个数为：", avg_N]
    f.writelines(list_cation)


if __name__ == "__main__":

    isInConfined_Cation(data)
    isInConfined_Anion(data)

    f.close()

效果如下：