prior_py.py 3.23 KB
Newer Older
1
2
3
4
5
import numpy as np
from scipy.spatial.distance import euclidean
from itertools import product
from ase.ga.utilities import closest_distances_generator

Malthe Kjær Bisbo's avatar
--    
Malthe Kjær Bisbo committed
6
class repulsive_prior():
7

Malthe Kjær Bisbo's avatar
--    
Malthe Kjær Bisbo committed
8
    def __init__(self, rcut=5, ratio_of_covalent_radii=0.7):
9
        self.rcut = rcut
Malthe Kjær Bisbo's avatar
--    
Malthe Kjær Bisbo committed
10
11
        self.ratio_of_covalent_radii = ratio_of_covalent_radii
        self.initialized = False
12

Malthe Kjær Bisbo's avatar
--    
Malthe Kjær Bisbo committed
13
    def initialize_from_atoms(self, atoms):
14
15
16
17
18
19
20
        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()
Malthe Kjær Bisbo's avatar
--    
Malthe Kjær Bisbo committed
21
22
23
24

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

        self.initialized = True
25
26

    def energy(self, a):
Malthe Kjær Bisbo's avatar
--    
Malthe Kjær Bisbo committed
27
28
29
        if not self.initialized:
            self.initialize_from_atoms(a)

30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
        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):
Malthe Kjær Bisbo's avatar
--    
Malthe Kjær Bisbo committed
48
49
50
        if not self.initialized:
            self.initialize_from_atoms(a)

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
80
81
82
83
84
85
86
87
88
89
90
91
        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