Create Bonds and Angles for LAMMPS's input

Posted 2020-09-09 花花今天没吃药

kkk

 1 #!/usr/bin/python
 2 #python 2.7
 3 text_all = []
 4 xatom_tmp = []
 5 xatom = []
 6 bond = []
 7 double_bond = []
 8 i = 0
 9 j = 0
10 dd = 0.0
11 #lattice constant
12 a = 10.86
13 b = 10.86
14 c = 10.86
15 
16 def isBond(list1, list2):
17 #i = -1 0 1,  j = -1 0 1
18     for i in range(-1,2):
19         for j in range(-1,2):
20             for k in range(-1,2):
21                 dd_tmp = 0.0
22                 dd_tmp += (list1[0] - list2[0] + a * i) ** 2
23                 dd_tmp += (list1[1] - list2[1] + b * j) ** 2
24                 dd_tmp += (list1[2] - list2[2] + c * k) ** 2
25                 dd_tmp = dd_tmp ** 0.5
26                 if dd_tmp < 3.0:
27                     return 1
28 
29     return 0
30 
31 def contain_same_point(list1, list2):
32     list3 = [0, 0.0, 0.0, 0.0]
33     for i in range(0,2):
34         for j in range(0,2):            
35             if(list1[i] == list2[j]):
36                 list3[0] = 1
37                 list3[1] = list1[1-i]
38                 list3[2] = list1[i]
39                 list3[3] = list2[1-j]
40     return list3
41                 
42 
43 file_in = open(‘data.lammps‘)
44 text_all = file_in.readlines()
45 file_in.close()
46 
47 num_lines = len(text_all)
48 #(0:5) means a[0], a[1], a[2], a[3], a[4], a[5], totally 5 elements
49 for i in range(19,83):
50     xatom_tmp.append(text_all[i].split())
51     xatom.append([xatom_tmp[i-19][2], xatom_tmp[i-19][3], xatom_tmp[i-19][4]])
52 
53 #initialize xatom_f as a 5x3 float matrix
54 xatom_f = [[0.0 for col in range(3)] for row in range(num_lines)] 
55 
56 for i in range(0,64):
57     xatom_f[i][0] = float(xatom[i][0])
58     xatom_f[i][1] = float(xatom[i][1])
59     xatom_f[i][2] = float(xatom[i][2])
60 
61 #bonds, A-B and B-A, same bond, count 1.
62 for i in range(0,64):
63     for j in range(i+1,64):
64         if (isBond(xatom_f[i], xatom_f[j]) == 1):
65             bond.append([i,j])
66 k = 0
67 print "Bonds: "
68 print "\n"
69 for i in bond:
70     k += 1
71     print k, 1, i[0], i[1]
72 
73 top = bond[0][0]
74 angle = []
75 print "Angles: "
76 print "\n"
77 for i in range(0,len(bond)):
78     for j in range(i+1,len(bond)):
79         tmp = contain_same_point(bond[i],bond[j])
80         if (tmp[0] == 1):
81             print tmp[1],tmp[2],tmp[3]