lammps统计六元环（非苯环）个数--Python实现

思路

1、六元碳环中，两原子最远距离为3X1.7=5.1

2、六个碳原子的集合中，每个碳原子彼此之间都只成两个C-C键的情况。只有一种可能————碳原子之间首尾相连连成六元环

步骤

1、找出距离目标原子距离<=5的所有原子，建立一个包含n个原子列表

2、在列表中随机取六个原子，建立一个小的合集，穷尽这种合集（排列组合C6/n）

3、哪一个合集中的六个原子满足彼此之间都恰好成两个键，则这一个合集能够组成一个苯环

#统计碳烟内部六元环的个数

import matplotlib.pyplot as plt
import numpy as np
import linecache
import itertools as it         #用于调用排列组合
from collections import Counter #用于调用列表中元素计数

#距离计算公式
def distant(l1,l2):
    a=np.sqrt((float(l1[2])-float(l2[2]))**2+(float(l1[3])-float(l2[3]))**2+(float(l1[4])-float(l2[4]))**2)
    return a 

def count_occurrences(lst, element):
    return Counter(lst)[element]

#打开碳烟颗粒的文件
with open("b.txt") as file_project:
    
    soot=[]
    lines =file_project.readlines()
    #print(lines)
    for line in lines:
        row = line.split()
        #print(row)
        soot.append(row)               #读取碳烟颗粒中每个原子信息
#print(soot)


C6=0
for atom1 in soot:
    circle=[]
    #circle.append(atom1)
    for atom2 in soot:
        if distant(atom1,atom2) <= 5:
            circle.append(atom2)
    #print(circle)
    combin=list(it.combinations(circle,6))
    print(len(combin))                  #排列组合的种类
    for num in combin:
        ever = []
        for atom1 in num:
            number=0
            for atom2 in num:
                if distant(atom1,atom2) <= 1.7:
                    number = number + 1
                    ever.append(number)                 
        a = count_occurrences(ever,2)
        print(a)
        if int(a) == 6 :
            C6 =C6 +1
print(C6)