-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathbrute.py
52 lines (41 loc) · 1.94 KB
/
brute.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
# Just documenting the problem setup.
print("N Lennard-Jones atoms are randomly places in a box.\n26 periodic images of each atom in this box are then generated.\nAverage distance and angle are calculated.\nFractional coordinate are used for distances and degress for angles.\nHopefully, this code is not meaningless!")
import numpy as n
N = int(input("Number of atoms you want to simulate: "))
# Generating random coordinates
r = n.random.rand(N, 3)
# A function that creates periodic images of each set of x-y-z coordinates
def PeriodicImages(R):
a = []
s = [-1, 0, 1]
for i in s:
x = R[0] + i
for j in s:
y = R[1] + j
for k in s:
z = R[2] + k
a.append([x, y, z])
coordinates = n.array(a)
return coordinates
# Testing the function is working as intended
atom = PeriodicImages(r[0])
print(atom)
print("Shape:", n.shape(atom))
# Here, and before calculating distances and angles, pair potential shall be calculated and minimized. I have been trying to do this for two weeks but the algorithms are just too long (and Allen & Tildesley are awful writers).
# Putting all coordinates in one list in preparation of distance and angle calculations
AllCoordinates = []
for i in range(n.shape(r)[0]):
AllCoordinates.extend(PeriodicImages(r[i]))
Distances = []
Angles = []
for i in AllCoordinates:
for j in AllCoordinates:
if (i != j).all(): # Doesn't calculate distances and angles for identical coordinates
d = n.sqrt( (i[0]-j[0])**2 + (i[1]-j[1])**2 + (i[2]-j[2])**2 )
Distances.append(d)
ang = n.degrees(n.arccos((n.dot(i, n.transpose(j))) / (n.linalg.norm(i)*n.linalg.norm(j)))) # Messy but correct
Angles.append(ang)
AverageDistance = sum(Distances)/len(Distances)
print("Average distance in fractional units =", AverageDistance)
AverageAngle = sum(Angles)/len(Angles)
print("Average angle in degrees =", AverageAngle)