prior_py.py 2.91 KB
Newer Older
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
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
import numpy as np
from scipy.spatial.distance import euclidean
from itertools import product
from ase.ga.utilities import closest_distances_generator

class delta():

    def __init__(self, atoms, rcut=5, ratio_of_covalent_radii=0.7):
        self.rcut = rcut

        self.pbc = atoms.get_pbc()
        self.cell = atoms.get_cell()

        self.cell_displacements = self.__get_neighbour_cells_displacement(self.pbc, self.cell)
        self.Ncells = len(self.cell_displacements)

        num = atoms.get_atomic_numbers()
        self.blmin = closest_distances_generator(num, ratio_of_covalent_radii=ratio_of_covalent_radii)

    def energy(self, a):
        num = a.get_atomic_numbers()
        x = a.get_positions()
        E = 0
        for i, xi in enumerate(x):
            for cell_index in range(self.Ncells):
                displacement = self.cell_displacements[cell_index]
                for j, xj in enumerate(x):
                    key = (num[i], num[j])
                    rmin = self.blmin[key]
                    radd = 1 - rmin
                    r = euclidean(xi, xj + displacement)
                    if r > 1e-5 and r < self.rcut:
                        r_scaled = r + radd
                        E += 1/r_scaled**12
        # Devide by two decause every pair is counted twice
        return E/2

    def forces(self, a):
        num = a.get_atomic_numbers()
        x = a.get_positions()
        Natoms, dim = x.shape
        dE = np.zeros((Natoms, dim))
        for i, xi in enumerate(x):
            for cell_index in range(self.Ncells):
                displacement = self.cell_displacements[cell_index]
                for j, xj in enumerate(x):
                    key = (num[i], num[j])
                    rmin = self.blmin[key]
                    radd = 1 - rmin
                    xj_pbc = xj+displacement
                    r = euclidean(xi,xj_pbc)
                    if r > 1e-5 and r < self.rcut:
                        r_scaled = r + radd
                        rijVec = xi-xj_pbc

                        dE[i] += 12*rijVec*(-1 / (r_scaled**13*r))
                        dE[j] += -12*rijVec*(-1 / (r_scaled**13*r))

        # Devide by two decause every pair is counted twice
        return - dE.reshape(-1) / 2

    def __get_neighbour_cells_displacement(self, pbc, cell):
        # Calculate neighbour cells
        cell_vec_norms = np.linalg.norm(cell, axis=0)
        neighbours = []
        for i in range(3):
            if pbc[i]:
                ncellmax = int(np.ceil(abs(self.rcut/cell_vec_norms[i]))) + 1  # +1 because atoms can be outside unitcell.
                neighbours.append(range(-ncellmax,ncellmax+1))
            else:
                neighbours.append([0])

        neighbourcells_disp = []
        for x,y,z in product(*neighbours):
            xyz = (x,y,z)
            displacement = np.dot(xyz, cell)
            neighbourcells_disp.append(list(displacement))

        return neighbourcells_disp